# Global and local stability for ghosts coupled to positive energy degrees of freedom

Cédric Deffayet, Aaron Held, Shinji Mukohyama, Alexander Vikman

► To cite this version:

Cédric Deffayet, Aaron Held, Shinji Mukohyama, Alexander Vikman. Global and local stability for ghosts coupled to positive energy degrees of freedom. Journal of Cosmology and Astroparticle Physics, 2023, 2023 (11), pp.031. 10.1088/1475-7516/2023/11/031. hal-04310024

HAL Id: hal-04310024

https://hal.science/hal-04310024v1

Submitted on 1 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
PAPER • OPEN ACCESS

# Global and local stability for ghosts coupled to positive energy degrees of freedom

To cite this article: Cédric Deffayet et al JCAP11(2023)031

View the [article online](#) for updates and enhancements.

You may also like

- A SO(1,3) gauge theory of quantum gravity: quantization and renormalizability proof

C Wiesendanger

- Lifshitz scaling, ringing black holes, and superradiance

Naritaka Oshita, Niayesh Afshordi and Shinji Mukohyama

- Breaking the warp barrier: hyper-fast solitons in Einstein-Maxwell-plasma theory

Erik W Lentz
<figure><img src="image_2.png" /></figure>

JCAP ANNIVERSARY
SPECIAL ISSUE

# Global and local stability for ghosts coupled to positive energy degrees of freedom

Cédric Deffayet<sup>D,a</sup> Aaron Held<sup>D,b,c</sup> Shinji Mukohyama<sup>d,e</sup> and Alexander Vikman<sup>D,f</sup>

<sup>a</sup>Laboratoire de Physique de l'École normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France

<sup>b</sup>Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universität Jena, Max-Wien-Platz 1, 07743 Jena, Germany

<sup>c</sup>The Princeton Gravity Initiative, Jadwin Hall, Princeton University, Princeton, New Jersey 08544, U.S.A.

<sup>d</sup>Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, 606-8502 Kyoto, Japan

<sup>e</sup>Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, 277-8583 Chiba, Japan

<sup>f</sup>CEICO — Central European Institute for Cosmology and Fundamental Physics, FZU — Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, 18221 Prague 8, Czech Republic

E-mail: cedric.deffayet@ens.fr, aaron.held@uni-jena.de, shinji.mukohyama@yukawa.kyoto-u.ac.jp, vikman@fzu.cz

**ABSTRACT:** Negative kinetic energies correspond to ghost degrees of freedom, which are potentially of relevance for cosmology, quantum gravity, and high energy physics. We present a novel wide class of stable mechanical systems where a positive energy degree of freedom interacts with a ghost. These theories have Hamiltonians unbounded from above and from below, are integrable, and contain free functions. We show analytically that their classical motion is bounded for all initial data. Moreover, we derive conditions allowing for Lyapunov stable equilibrium points. A subclass of these stable systems has simple polynomial potentials with stable equilibrium points entirely due to interactions with the ghost. All these findings are fully supported by numerical computations which we also use to gather evidence for stability in various nonintegrable systems.

**KEYWORDS:** dark energy theory, dark matter theory

ARXIV ePRINT: 2305.09631

<figure><img src="image_3.png" /></figure>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Integrable interacting ghosts via complex canonical transformations</td><td>3</td></tr><tr><td>3</td><td>The integrable Liouville Model &amp; a proof of bounded motion</td><td>7</td></tr><tr><td>3.1</td><td>Proof of boundedness of motion</td><td>9</td></tr><tr><td>3.2</td><td>A first application of the proof</td><td>14</td></tr><tr><td>3.3</td><td>A polynomial subclass of the Liouville model: towards more general conditions for stable motion</td><td>15</td></tr><tr><td>3.4</td><td>An unstable polynomial subclass: refining more general conditions for stable motion</td><td>20</td></tr><tr><td>3.5</td><td>Lagrange vs Lyapunov stability</td><td>23</td></tr><tr><td>4</td><td>Numerical investigation of nonintegrable models</td><td>30</td></tr><tr><td>4.1</td><td>Numerical setup</td><td>31</td></tr><tr><td>4.2</td><td>Benchmark: the Liouville model revisited</td><td>32</td></tr><tr><td>4.3</td><td>Numerical evidence in the absence of integrability</td><td>34</td></tr><tr><td>5</td><td>Discussion</td><td>37</td></tr><tr><td>A</td><td>Integrable two-particle Hamiltonians with polynomial constants of motion up to quadratic order in the momenta</td><td>39</td></tr><tr><td>B</td><td>Integrability and vacua</td><td>41</td></tr><tr><td>C</td><td>Stable motion with a ghostly x</td><td>43</td></tr><tr><td>D</td><td>Polynomial Liouville potentials</td><td>45</td></tr><tr><td>E</td><td>Polynomial class-3 potentials</td><td>47</td></tr></table>

## 1 Introduction

Ghosts are dynamical degrees of freedom (DoF), described by a Hamiltonian with negative kinetic energies, as opposed to the kinetic energies of standard degrees of freedom. Such ghosts generically appear in the context of theories with higher derivatives [1] (see also [2–4] for reviews) and have been considered in a variety of physical scenarios: they have been advocated as a way to regulate the ultraviolet behaviour of quantum field theories [5–8]; they have been suggested as a way to address the cosmological constant problem [9–11]; they have been used to obtain bouncing cosmologies [12]; and they have been proposed in the context of dark energy [13] or, more recently, as a way to ameliorate the Hubble tension [14].
In spite of these applications, ghosts are typically disregarded since they are expected to generically lead to catastrophic instabilities, both at the classical and at the quantum level (see e.g [15]) — for the purpose of this work, we will only be concerned with the classical case.

Whenever positive- and negative-kinetic-energy degrees of freedom do not interact, the presence of a ghost is trivially harmless. Without coupling, no energy exchange can take place and the respective ghostly and non-ghostly Hamiltonians can each be stable due to suitable self-interactions.

In contrast, it is widely expected that a generic coupling between the positive- and negative-kinetic-energy degrees of freedom, leads to an energy exchange between them that instigates an unbounded growth of absolute values of both energies. Hence, one expects a runaway evolution to larger and larger phase-space variables. However, in previous work [16], coauthored by some of the authors of the current paper, a specific counterexample to such a catastrophe was presented. Besides completing a stability proof for a class of models introduced in [16], the key physical motivation for this work can be understood in terms of the following simple question: is the stable model from [16] too specific or can stability be extended to a wider class of integrable and nonintegrable systems with ghosts?

The potential for stable motion in the presence of interacting ghosts has, of course, been previously considered. (While we will be concerned only with classical theories, we refer to [6, 8, 17–27] for proposals at the quantum level.) At the classical level, it has been argued that theories with ghosts can be considered “benign”, provided that the runaway is sufficiently mild and the motion can be extended to the infinite future [25, 28–31]. At the same time, indications have been found that classical motion in the presence of ghosts can exhibit so-called “islands of stability” (i.e., bounded motion for a restricted set of initial conditions) [21, 28–36]. All of these indications are numerical (and, as such, cannot be fully conclusive, as they cannot cover the whole evolution of the system to the infinite future along even a single trajectory, not to mention all the Hamiltonian trajectories) and/or hold only for restricted sets of initial conditions, i.e., do not conclusively advocate for global stability. Moreover, we still lack a physical rationale or criterion to distinguish potentially stable ghostly interactions from unstable ones.

As mentioned before, the recent work [16] has established the first conclusive proof of globally stable motion for a positive-energy harmonic oscillator interacting through a specific non-polynomial potential with a ghost in the form of a negative-energy oscillator. It was shown that, for arbitrary initial conditions, the motion remains bounded. The proof relies on the integrable nature of the model. We stress that this kind of global stability (also referred to as the Lagrange stability) is different from the above notions of “benign” ghosts and/or “islands of stability”. We also stress that there is no direct connection of global stability to local (Lyapunov) stability [37]. In this context, we would like to mention that [16] has also proven Lyapunov stability of the equilibrium configuration of the model. It is global (Lagrange) stability that we are mostly concerned with and we will, from here on, refer to it simply as stability, everywhere where this would not cause confusion. One of the reasons for this concern is that quantum mechanics effectively probes the whole configuration space. Thus, regions beyond the, usually small, “islands of stability” can force the wave function and probability distributions to runaway to large values of phase-space coordinates in the quantized theory.
The purpose of this paper is the following: first, we extend the proof in [16] to a much larger class of integrable models with kinetic terms quadratic in canonical momenta and with two free functions of (particular combinations) of the coordinates canonically conjugated to them. This class of systems encompasses, in particular, a tower of stable polynomial interactions. Our proof also covers the subclass that has been announced in [16]. Second, we prove that a wide subclass of these models possesses locally stable vacua. Third, we propose physical criteria for stable motion of more generic systems with ghosts interacting with positive energy degrees of freedom. Fourth, we investigate the above numerically, both for integrable and nonintegrable models.

With these goals in mind, we start in section 2 by introducing a generic procedure to relate classical mechanical systems of $N$ degrees of freedom with and without ghosts by means of suitable complex canonical transformations. We then consider a particular class of integrable models in section 3, for which in section 3.1 we provide a fully analytic proof of the boundedness of motion, identify a polynomial subclass in section 3.3, and discuss Lyapunov stability of equilibrium points in section 3.5. In particular, we find that new Lyapunov stable vacua can be instigated by polynomial interactions with a ghost. Incidentally, we stress that the polynomial models introduced in section 3.3 stay polynomial and integrable when the degrees of freedom all have positive kinetic energies and appear as such to be novel to this work. Finally, we extend the discussion to nonintegrable models in section 4, for which we gather numerical evidence that stable (or at least longlived) motion can persist in absence of integrability. We close with our conclusions in section 5. We provide additional technical material in appendices. In particular, in appendix A we list various classes of integrable models; in appendix B we derive equations defining the models with the first integral quadratic in canonical momenta and demonstrate that Lyapunov stable equilibrium points are located in the saddle points of the potential; in appendix C we discuss stability in alternatively complexified Liouville models; in appendix D and appendix E we provide more details on integrable models with polynomial potentials.

## 2 Integrable interacting ghosts via complex canonical transformations

We start by considering $N$ classical mechanical degrees of freedom (DoFs) with a Hamiltonian

$$H(x_i, p_i) = \sum_{i=1}^{N} \frac{p_i^2}{2m_i} + V(x_i, \dots, x_N), \qquad (2.1)$$

where $p_i$ denotes the canonical momenta conjugate to the coordinates $x_i$. Together, these denote the phase-space variables $\xi_i = (x_i, p_i)$. For now, all $N$ DoFs have same-sign kinetic terms and thus correspond to modes with positive kinetic energy, as the masses $m_i$ are assumed to be positive. We notice that a rotation in the complex plane to the imaginary axis of any of the momentum variables, i.e., $p_n \to \pm ip_n$, flips the sign of the respective kinetic term, i.e., $p_n^2/m_n \to -p_n^2/m_n$, and thus turns the n-th degree of freedom (DoF) with positive kinetic energy into a DoF with negative kinetic energy. Hence, this transformation effectively flips the sign of the mass $m_n$. This rotation in the complex plane can be completed into a complex canonical transformation by simultaneously transforming the respective $x_n \to \mp ix_n$.
The resulting (univalent) transformation

$$ \mathcal{C}_{x^n}^{\pm} : \left\{ \begin{array}{l} x_j \rightarrow \pm i x_j \quad \text{for } j = n \\ p_j \rightarrow \mp i p_j \quad \text{for } j = n \\ x_j \rightarrow x_j \quad \forall j \neq n \\ p_j \rightarrow p_j \quad \forall j \neq n \end{array} \right\}, \qquad (2.2) $$

is canonical in the sense that it preserves the Poisson bracket $\{\}$, i.e.,

$$ \begin{align} \{\mathcal{C}_{x^n}^{\pm}(x_i), \mathcal{C}_{x^n}^{\pm}(p_j)\}_{\text{PB}} &= -\{\mathcal{C}_{x^n}^{\pm}(p_i), \mathcal{C}_{x^n}^{\pm}(x_j)\}_{\text{PB}} = \delta_{ij}, \nonumber \\ \{\mathcal{C}_{x^n}^{\pm}(x_i), \mathcal{C}_{x^n}^{\pm}(x_j)\}_{\text{PB}} &= 0, \nonumber \\ \{\mathcal{C}_{x^n}^{\pm}(p_i), \mathcal{C}_{x^n}^{\pm}(p_j)\}_{\text{PB}} &= 0. \tag{2.3} \end{align} $$

Given that $\mathcal{C}_{x^n}^{\pm}$ is a canonical transformation, it moreover preserves the Hamilton equations, the time-independence of the Hamiltonian itself, as well as all the time-independence of any other existing constant of motion. We note also that $\mathcal{C}_{x^n}^{\pm} \circ \mathcal{C}_{x^n}^{\mp} = \mathcal{C}_{x^n}^{\mp} \circ \mathcal{C}_{x^n}^{\pm} = 1$. Moreover, any function $F(\xi)$ of the phase-space variables $\xi$ which is even under the simultaneous change of the sign of the variables $x_n$ and $p_n$, i.e.,

$$ F(x_{i \neq n}, x_n, p_{i \neq n}, p_n) = F(x_{i \neq n}, -x_n, p_{i \neq n}, -p_n), \qquad (2.4) $$

is invariant under the successive application of two identical $\mathcal{C}_{x^n}^{\pm}$. We will also refer to such functions as ‘$(x_n, p_n)$-parity even’. On the other hand, we refer to a potential satisfying $V(x_{i \neq n}, x_n) = V(x_{i \neq n}, -x_n)$ as ‘$x_n$-parity even’. Obviously, the Hamiltonian in eq. (2.1) is $(x_n, p_n)$-parity even if and only if the potential is $x_n$-parity even.

In general, $\mathcal{C}_{x^n}^{\pm}$ does not preserve real-valuedness of the Hamiltonian. (The same, and all of what follows, also holds for other constants of motion.) However, real-valuedness is obviously preserved if the potential is polynomial (or Taylor expandable with respect to $x_n$) and $x_n$-parity even. This is the case, e.g., for the canonical Hamiltonian of $N$ noninteracting positive kinetic-energy modes (where here and henceforth we normalize the kinetic terms just to be $\pm p_i^2/2$, i.e., we set the masses to one, but keep the possibility to have a different frequency $\omega_i$ for each mode)

$$ H_{\text{canonical}} = \sum_{i=1}^{N} \frac{1}{2} (p_i^2 + \omega_i^2 x_i^2). \qquad (2.5) $$

In this case, the canonical transformation $\mathcal{C}_{x^n}^{\pm}$ transforms a real-valued Hamiltonian with $N$ positive-energy modes to a dual real-valued Hamiltonian with $N-1$ positive kinetic-energy DoFs and one negative kinetic-energy DoF, i.e., a ghost. This can easily be extended by performing several distinct such canonical transformations to introduce multiple ghost degrees of freedom. For the general Hamiltonian in eq. (2.1), if the potential is parity-even in all of the respective variables, then all of these Hamiltonians are real-valued as far as the original one is. All of the above also holds if further parity-even interactions are included.

In the rest of this work, we will specify to systems with two degrees of freedom $x$ and $y$ and with a Hamiltonian

$$ H = \frac{1}{2} p_x^2 + \sigma \frac{1}{2} p_y^2 + V(x, y), \qquad (2.6) $$
where $\sigma = \pm 1$ decides between the case of two positive kinetic terms ($\sigma = +1$) and the case of two opposite-sign kinetic terms ($\sigma = -1$). We refer to $\sigma = +1$ as the “PP” case, denoting two positive kinetic-energy DoFs, and to $\sigma = -1$ as the “PG” case, where a positive kinetic-energy DoF is coupled to a negative kinetic-energy (ghost) DoF. Starting with a PP model, and applying a complex canonical transformation (cf. eq. (2.2)) on the variable $y$, we obtain a new model of PG type with Hamiltonian $\frac{1}{2}p_x^2 - \frac{1}{2}p_y^2 + V(x, \pm iy)$ which can be real depending on the $V(x, y)$ one starts out with.

In the following, it will be sometimes convenient to split the potential $V(x, y)$ into three parts as follows for a PG system with $\sigma = -1$,

$$V_P(x) = V(x, 0), \quad \text{and} \quad V_G(y) = V(0, y), \quad \text{while} \\ V_{\text{int}}(x, y) = V(x, y) - V_P(x) - V_G(y), \qquad (2.7)$$

such that $V = V_P + V_G + V_{\text{int}}$. We refer to $V_P(x)$ and $V_G(y)$ as the decoupled potentials, to $V_{\text{int}}(x, y)$ as the interaction potential, and to $V(x, y)$ as the full potential. This makes sense, in particular, if the structure of the potential $V(x, y)$ is such that $y = 0$ is a solution of $y$-equation of motion for arbitrary $x(t)$ and analogously $x = 0$ is a solution of the $x$-equation of motion for arbitrary $y(t)$. If this property holds,<fnref n="1" /> one can separate self-interactions from cross-coupling interactions.

More specifically, for what concerns the analytic part of this work, our starting point will be integrable PP systems with a time independent Hamiltonian $H$ of the form in eq. (2.6) and with an extra constant of motion $I(\xi)$ such that

$$\frac{dI}{dt} = \{I, H\}_{\text{PB}} = 0. \qquad (2.8)$$

The study of such systems has been pioneered by Liouville [38, 39] who showed that if an $N$ dimensional Hamiltonian system possesses $N$ functionally independent and Poisson-commuting constants of motion, then the motion is integrable (by quadrature). In principle, there can be more than $N$ (i.e., up to $(2N-1)$) functionally independent constants of motion. If so, the system is referred to as superintegrable.

Some Hamiltonian systems are trivially integrable. For instance, the canonical Hamiltonian of $N$ uncoupled harmonic oscillators, see eq. (2.5), possesses $N$ conserved quantities: the individual Hamiltonians of each oscillator. For other Hamiltonian systems, integrability is, of course, nontrivial and takes more effort to be revealed.

One way to systematically search for constants of motion is by direct methods. Direct methods (i) make an ansatz for $I(\xi)$ — typically a polynomial in momenta, dressed with undetermined functions as prefactors — and then (ii) classify solutions to the resulting (highly overdetermined) set of partial differential equations (PDEs) implied by $\{I, H\}_{\text{PB}} = 0$. When restricting to two particle systems, such direct searches for constants of motion (and the corresponding integrable potentials) have been performed at fixed polynomial order in the momenta. Up to quadratic order, the classification has been pioneered in [40], formalized

<footnote n="1">¹It may be that this property does not hold for the pair $(x, y)$ but holds for a pair $(x', y')$ that is related to $(x, y)$ by a Poincaré transformation. Then, one can redefine variables as $(x, y) \rightarrow (x', y')$. Of course, this property may not hold at all.</footnote>
in [41], and subsequently completed in [42–49]. We summarize the complete set of resulting pairs of integrable potentials $V(x, y)$ and constants of motion in appendix A. At cubic order in the momenta and beyond, no full classification is known, cf. [49] for a comprehensive review of partial results and for references to related methods.

As stressed above, the on-shell constancy of some quantity $I$ is preserved by the action of $\mathcal{L}_{x^n}$ as a mere consequence of the algebraic nature of equation eq. (2.8) and the canonical nature of the transformation $\mathcal{L}_{x^n}$. Hence, acting on a PP integrable system with a transformation $\mathcal{L}_{x^n}$, we obtain a PG integrable one. Upon this action, the Hamiltonian can, in general, become complex. However, an inspection of the potentials given in appendix A shows that there exist large classes of PP integrable systems yielding real-valued PG integrable systems upon application of the complex canonical transformation $\mathcal{L}_y$ (or $\mathcal{L}_x$). This is the case, e.g., for class 1 of appendix A. In the next section, we use this real-valued PP $\Leftrightarrow$ PG pair to analytically discuss the conditions for bounded motion in the presence of a ghost.

In passing, we also note that there can be real-valued integrable PG systems, for which the PP side is complex-valued, and which may thus have not yet received much attention. An integrable example of this kind is given by class 5 in appendix A. In this case, the application of the complex canonical transformation $\mathcal{L}_y$ (or $\mathcal{L}_x$) transforms a complex $V$ to a real one.

To avoid confusion, we stress that complex canonical transformations $\mathcal{L}_{x^n}$ preserve the physics only when accompanied by the corresponding rotations in the complex planes of the initial data. Otherwise, the transformed theory is physically different and equivalent to probing the original theory with purely imaginary initial data. Therefore, the dynamical properties of such “ghostified” theories can be completely different from those of the original theory when all quantities considered are real, which we will assume in our PG models.

Let us illustrate this with a simple example of a theory with one real degree of freedom $x$ (and real $p$)

$$H = \frac{p^2}{2} + \frac{1}{4} \cosh x. \qquad (2.9)$$

Clearly, for arbitrary, we stress real, initial data, the solution perpetually evolves in a finite region of phase space bounded by the potential well $\cosh x$ and the positive kinetic energy. This solution is a closed trajectory in the phase space $(x, p)$, see the upper panel of figure 1. Now let us perform the canonical transformation $\mathcal{L}_x$, so that $p = i\bar{p}$ and $x = -i\bar{x}$, transforming eq. (2.9) to the new theory of real variables $(\bar{x}, \bar{p})$ and Hamiltonian, i.e.,

$$\bar{H} = -\frac{\bar{p}^2}{2} + \frac{1}{4} \cos \bar{x}. \qquad (2.10)$$

We could consider the solution $(\bar{x}(t), \bar{p}(t))$ for the system in eq. (2.10), with purely imaginary initial data $(\bar{x}_0, \bar{p}_0) = (ix_0, -ip_0)$ to obtain just the real solution $(x(t), p(t)) = (-i\bar{x}(t), i\bar{p}(t))$ of eq. (2.9) with the real initial data $(x_0, p_0)$. In particular, $H = \bar{H}$ for all the solutions which are related in this way. However, a solution of eq. (2.10) with real initial data $(\bar{x}_0, \bar{p}_0)$ corresponds to a purely imaginary solution of eq. (2.9) with initial data $(-i\bar{x}_0, i\bar{p}_0)$. Thus, purely real solutions of eq. (2.9) are very different from purely real solutions of eq. (2.10) even though both systems are related by a simple canonical transformation. In particular, for arbitrary real initial coordinates $\bar{x}_0$ and initial momenta $\bar{p}_0 > 1$ (or $\bar{p}_0 < -1$), the trajectory
<figure><img src="image_4.png" /><figcaption>**Figure 1.** We show the phase space $(x, p)$ of the original Hamiltonian in eq. (2.9) (upper panel) as well as the phase space $(\bar{x}, \bar{p})$ of the transformed Hamiltonian in eq. (2.10) (lower panel). In the ghostified case (lower panel), all trajectories with initial momentum $|\bar{p}_0| > 1$ (cf. red horizontal lines) evolve over the potential barrier and run away to plus or minus infinity in $\bar{x}$, irrespective of the initial value of $\bar{x}_0$. Two such runaway trajectories are depicted in blue. In contrast, all phase space trajectories of eq. (2.9) (upper panel) are closed curves.</figcaption></figure>

of eq. (2.10) evolves over the potential barrier and is unbounded, with a runaway of $\bar{x}(t)$, see the right panel of figure 1. As we mentioned earlier, there are no such real-valued trajectories in eq. (2.9). Clearly, in the phase space $(\bar{x}, \bar{p})$ of eq. (2.10) such real-valued trajectories are not closed, whereas all real-valued trajectories in $(x, p)$ of eq. (2.9) are closed, thus even the topology of the phase space trajectories is changed. In this regard it is worth mentioning that with real $(\bar{x}, \bar{p})$ there are not only an infinite number of new stable equilibrium points, but also new unstable equilibrium points in-between the former, see figure 1. While the above refers to local (Lyapunov) stability, the “ghostified” real theory in eq. (2.9) is globally (Lagrange) unstable, in contrast to the original real theory eq. (2.9). Two of the respective runaway trajectories are marked as thick blue lines in figure 1.

### 3 The integrable Liouville Model & a proof of bounded motion

In the following, we consider an integrable class of theories, cf. class 1 in appendix A, for which we can rigorously establish boundedness of motion. Specifically, we focus on

$$H_{LV} = \frac{p_x^2}{2} + \sigma \frac{p_y^2}{2} + V_{LV}(x, y) \qquad (3.1)$$

$$\text{with} \quad V_{LV} = \frac{f(u) - g(v)}{u^2 - v^2}, \qquad (3.2)$$
where $f$ and $g$ are arbitrary functions of the coordinates $(u, v)$. The coordinates $(u, v)$ are related to $(x, y)$ via

$$u^2 = \frac{1}{2} \left( r^2 + c + \sqrt{(r^2 + c)^2 - 4cx^2} \right), \qquad (3.3)$$

$$v^2 = \frac{1}{2} \left( r^2 + c - \sqrt{(r^2 + c)^2 - 4cx^2} \right), \qquad (3.4)$$

$$r^2 = x^2 + \sigma y^2, \qquad (3.5)$$

where $c$ is an arbitrary constant. Besides the Hamiltonian, this theory exhibits a second constant of motion, i.e.,

$$I_{LV} = -\sigma(p_y x - \sigma p_x y)^2 - c p_x^2 + \mathcal{V} \qquad (3.6)$$

$$\text{with } \mathcal{V} = 2 \frac{u^2 g(v) - v^2 f(u)}{u^2 - v^2}, \qquad (3.7)$$

where we have (for later purpose) defined the momentum-independent part of the constant of motion $\mathcal{V} = I_{LV}|_{p_x, y \to 0}$. Although this model was initially studied in the PP case (i.e. with $\sigma = +1$) the discussion of the previous section shows that it stays integrable also for $\sigma = -1$, i.e., in the PG case. The corresponding extra integral of motion (i.e. eq. (3.6) with $\sigma = -1$) is obtained from the PP one (i.e. eq. (3.6) with $\sigma = +1$), as explained in the previous section.

For $\sigma = 1$ and $c \ge 0$, $u$ and $v$ define so-called elliptic coordinates, i.e., curves of constant $u \in (\sqrt{c}, \infty)$ trace out ellipses and curves of constant $v \in (0, \sqrt{c})$ trace out hyperbolas in the $(x, y)$ plane. For $\sigma = -1$ and $c \le 0$, curves of constant $u \in (0, \infty)$ and curves of constant $-iv \in (\sqrt{c}, \infty)$ both trace out hyperbolas in the $(x, y)$ plane, cf. figure 2.

We note that, even though this class of integrable systems was derived, and is usually presented, with an arbitrary constant $c$, it is only the sign of this constant which influences the dynamics. Indeed, for $c \ne 0$ one can first perform a canonical transformation

$$\begin{aligned} x &\rightarrow |c|^{1/2} x, & y &\rightarrow |c|^{1/2} y, \\ p_x &\rightarrow |c|^{-1/2} p_x, & p_y &\rightarrow |c|^{-1/2} p_y. \end{aligned} \qquad (3.8)$$

Then one can complete it with a rescaling of time

$$t \rightarrow |c|t, \qquad H \rightarrow |c|^{-1} H. \qquad (3.9)$$

These transformations completely absorb $|c|$ into the arbitrariness of $f(u)$ and $g(v)$, except for the limit cases $c \to \infty$ or $c \to 0$. For instance, the latter limit transforms the Liouville model to a subclass of class 2, see eq. (A.9), with $f = 0$.

The above integrable Hamiltonian in eq. (3.1) follows as a subclass of a system first investigated by Liouville [50] (in the PP case). Liouville starts out from a Hamiltonian with non-standard kinetic terms, i.e.,

$$H_{LV} = \frac{1}{F(\alpha) - G(\beta)} \left[ \frac{p_{\alpha}^2}{2} + \frac{p_{\beta}^2}{2} + (f(\alpha) - g(\beta)) \right], \qquad (3.10)$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left-hand panel: curves of constant $u$ (thick) and $v$ (thin) in the $(x, y)$ plane for the Liouville model without ghost ($\sigma = +1$). Right-hand panel: curves of constant $u$ (thick) and $-iv$ (thin) in the $(x, y)$ plane for the Liouville model with ghost ($\sigma = -1$).</figcaption></figure>

where $(\alpha, \beta)$ are yet another set of coordinates, $(p_\alpha, p_\beta)$ are the respective canonical momenta, and $f(\alpha)$, $F(\alpha)$, $g(\beta)$, $G(\beta)$ are arbitrary functions. Following Liouville, we found that

$$I_{LV} = p_{\alpha}^{2} + 2f(\alpha) - 2F(\alpha) H_{LV} \qquad (3.11)$$

Poisson-commutes with the Hamiltonian $H_{LV}$ and thus defines a second constant of motion.<fnref n="2" /> Hence, the model is integrable, irrespective of the choice of $f(\alpha)$, $F(\alpha)$, $g(\beta)$, and $G(\beta)$.

The coordinates $(\alpha, \beta)$ are related to $(u, v)$ via

$$\begin{aligned} u^{2} &= c \sinh \alpha^{2}, \\ v^{2} &= c \sin \beta^{2}. \end{aligned} \qquad (3.12)$$

In coordinates $(u, v)$, the Liouville Hamiltonian $H_{LV}$ and constant of motion $I_{LV}$ read

$$H_{LV} = \frac{1}{F-G} \left[ \frac{p_u^2}{2}(u^2 - c) + \frac{p_v^2}{2}(c - v^2) + (f-g) \right], \qquad (3.13)$$

$$I_{LV} = \frac{1}{F-G} \left[ \frac{p_u^2}{2}(u^2 - c)G + \frac{p_v^2}{2}(c - v^2)F - (fG - gF) \right]. \qquad (3.14)$$

For the specific choice of $F(u) = u^2$ and $G(v) = v^2$, and with the transformation from $(u, v)$ to $(x, y)$, these reduce to the Hamiltonian and constant of motion previously given in eqs. (3.1) and (3.6) with $\sigma = +1$. For more general $F(u)$ and $G(v)$, the transformation to $(x, y)$ coordinates does not bring the Hamiltonian to a form with standard kinetic terms. In the following, we will thus restrict to the case of $F(u) = u^2$ and $G(v) = v^2$.

## 3.1 Proof of boundedness of motion

We show here that the PG Liouville model with $\sigma = -1$ (in the form of eq. (3.1) and with the additional constant of motion in eq. (3.6)) is such that — under certain conditions on $f(u)$ and $g(v)$ detailed below — the phase-space motion is always bounded, irrespective of

<footnote n="2">²Equivalently, $K_{LV} = p_{\beta}^2 - 2g(\alpha) + 2G(\alpha) H_{LV}$ denotes a constant of motion, but $K_{LV} = -I_{LV}$ and is therefore not independent.</footnote>
the choice of initial conditions. This applies to the model considered in [16], which is a special class of PG Liouville models, but also extends the result of [16] to a much larger class of ghostly oscillator interacting with a positive energy oscillator where the motion can hence be proven analytically not to run away.

We start by defining three further helpful coordinate combinations

$$ 
\begin{aligned}
W &\equiv u^2 + v^2 = x^2 - y^2 + c \\
  &= r^2 + c,
\end{aligned} 
\qquad (3.15) 
$$

$$ 
\begin{aligned}
\tilde{W} &\equiv u^2 - v^2 \\
&= \sqrt{r^4 - 2c(x^2 + y^2) + c^2},
\end{aligned} 
\qquad (3.16) 
$$

$$ w \equiv -4 u^2 v^2 = -4 c x^2, \qquad (3.17) $$

with which we may re-express $\tilde{W} = \sqrt{W^2 + w}$ and write

$$ u^2 = \frac{1}{2} (W + \sqrt{W^2 + w}) = \frac{1}{2} (W + \tilde{W}) \qquad (3.18) $$

$$ v^2 = \frac{1}{2} (W - \sqrt{W^2 + w}) = \frac{1}{2} (W - \tilde{W}) \qquad (3.19) $$

We choose here and henceforth a negative $c < 0$ such that the momentum-dependent part of the constant of motion (cf. the first two terms in eq. (3.14)) are positive. Furthermore, a negative $c < 0$ implies

$$ v^2 \le 0, \qquad (3.20) $$

and while the above $v^2$ (given in eq. (3.18)) is a real expression, the so defined $v$ is purely imaginary. In order to manipulate real and positive variables, we define the real and positive quantities $\tilde{v}$ and $\tilde{c}$ as

$$ \tilde{v} = |v| \qquad (3.21) $$

$$ \tilde{c} = |c| = -c \qquad (3.22) $$

such that $v^2 = -\tilde{v}^2$, and thus we have

$$ u^2 \ge 0, \qquad (3.23) $$

$$ \tilde{v}^2 \ge 0, \qquad (3.24) $$

$$ \tilde{W} \ge 0, \qquad (3.25) $$

$$ \tilde{c} > 0, \qquad (3.26) $$

all these quantities being real, as well as, by choice, the functions $f$ and $g$.

We now prove that phase-space motion is bounded if

(i) $f(u)$ and $g(v)$ are bounded below, i.e.,

$$ f(u) \ge f_0, \qquad (3.27) $$

$$ g(v) \ge g_0, \qquad (3.28) $$

with constants $f_0, g_0 \in \mathbb{R}$; and
(ii) at large $|u|$ and $|v|$, these lower bounds sharpen to

$$f(u) \geq 4F_0 |u|^{\zeta} > 0, \qquad (3.29)$$

$$g(v) \geq 4G_0 |v|^{\eta} > 0, \qquad (3.30)$$

with positive constants $F_0, G_0 \in \mathbb{R}^+$ as well as $\zeta > 2$ and $\eta > 2$ (and the factor 4 is just here for later convenience).

The first step of our proof is the introduction of a new first integral $J_{LV}$ defined as (using eqs. (3.21), (3.24) and (3.26))

$$\begin{aligned} J_{LV} &= I_{LV} - \tilde{c} H_{LV} \\ &= (x p_y + y p_x)^2 + \frac{\tilde{c}}{2} (p_x^2 + p_y^2) + \mathcal{U}, \end{aligned} \qquad (3.31)$$

where

$$\begin{aligned} \mathcal{U}(u, v) &= \frac{(2 u^2 + \tilde{c})g(v) + (2 \tilde{v}^2 - \tilde{c})f(u)}{u^2 + \tilde{v}^2} \\ &= \gamma_+ g(v) + \gamma_- f(u), \end{aligned} \qquad (3.32)$$

and $\gamma_{\pm}$ are obtained to be (using eqs. (3.16), (3.18) and (3.19))

$$\gamma_{\pm} = 1 \pm \frac{r^2}{\bar{W}}. \qquad (3.33)$$

Using eq. (3.16), we get $-1 < \frac{r^2}{\bar{W}} < 1$ implying

$$0 < \gamma_{\pm} < 2. \qquad (3.34)$$

This, using the condition (i) above implies that

$$\gamma_+ g(v) \geq -2|g_0|, \qquad (3.35)$$

$$\gamma_- f(u) \geq -2|f_0|, \qquad (3.36)$$

where the equality may hold if $g_0 = 0$ or/and $f_0 = 0$, and hence, using eq. (3.32), we get that

$$\mathcal{U} \geq -2(|g_0| + |f_0|), \qquad (3.37)$$

showing that $\mathcal{U}$ is bounded below. As $J_{LV}$ defined in eq. (3.31) is conserved, this implies that $|x p_y + y p_x|, |p_x|$ and $|p_y|$ are bounded on shell. This is the first step of our proof. As it should be clear, this first step puts forward an analogy between a conserved and bounded from below Hamiltonian and the conserved $J_{LV}$, where all the “kinetic energies” are positive definite and the “potential energy” $\mathcal{U}(u, v)$ has just been shown to be bounded from below.

To proceed further, we use again this analogy, showing that conditions (ii) imply that $\mathcal{U}(x, y)$ is growing without bound at large radius $R = \sqrt{x^2 + y^2}$, which excludes runaways of the dynamical variables $x$ and $y$. To do so we introduce polar coordinates $(R, \phi)$ such
that $x = R \cos \phi$ and $y = R \sin \phi$, yielding $r^2 = x^2 - y^2 = R^2 \kappa$ where we have defined $\kappa = \cos 2\phi$. As a consequence, one has

$$ \gamma_{\pm} = 1 \pm \frac{R^2 \kappa}{\sqrt{R^4 \kappa^2 + 2R^2 \tilde{c} + \tilde{c}^2}}. \qquad (3.38) $$

As we are interested in large $R$, we can assume that

$$ \sqrt{2R} \gg \sqrt{\tilde{c}}, \qquad (3.39) $$

implying, in particular, that the last term in the square-root in eq. (3.38) is much smaller than the second one. We also have from eq. (3.16)

$$ u^2 + \tilde{v}^2 = \sqrt{R^4 \kappa^2 + 2R^2 \tilde{c} + \tilde{c}^2} > \sqrt{R^4 \kappa^2 + 2R^2 \tilde{c}}. \qquad (3.40) $$

Using this we conclude that one has either

$$ u^2 > \frac{1}{2} \left( \sqrt{R^4 \kappa^2 + 2R^2 \tilde{c}} \right), \qquad (3.41) $$

or

$$ \tilde{v}^2 > \frac{1}{2} \left( \sqrt{R^4 \kappa^2 + 2R^2 \tilde{c}} \right). \qquad (3.42) $$

Distinguishing the above, and moreover between

$$ R^4 \kappa^2 < 2R^2 \tilde{c} \Leftrightarrow R|\kappa| < \sqrt{2\tilde{c}}, \qquad (3.43) $$

$$ \text{or} \quad R^4 \kappa^2 \ge 2R^2 \tilde{c} \Leftrightarrow R|\kappa| \ge \sqrt{2\tilde{c}}, \qquad (3.44) $$

we will now show that, in all the respective four cases, $\mathcal{U}$ grows at least as fast as a suitable power of $R$.

Let us first assume that eq. (3.43) holds. In this case, we have

$$ \sqrt{R^4 \kappa^2 + 2R^2 \tilde{c} + \tilde{c}^2} > \sqrt{2R^4 \kappa^2}. \qquad (3.45) $$

This implies a lower bound on $\gamma_{\pm}$ given by<fnref n="3" />

$$ \gamma_{\pm} \ge 1 - \frac{1}{\sqrt{2}} > \frac{1}{4}, \qquad (3.46) $$

which holds for either sign of $\kappa$ (and including $\kappa = 0$). Moreover, we have from eqs. (3.41) and (3.42) that either $u^2 \ge R\sqrt{\tilde{c}/2}$ or $\tilde{v}^2 \ge R\sqrt{\tilde{c}/2}$. In the first subcase, we can use the bounds eqs. (3.29), (3.35) and (3.46) to get that at large $R$

$$ \mathcal{U} \ge -2|g_0| + F_0 (\tilde{c}/2)^{\zeta/4} R^{\zeta/2}. \qquad (3.47) $$

In the second subcase, we get similarly

$$ \mathcal{U} \ge -2|f_0| + G_0 (\tilde{c}/2)^{\eta/4} R^{\eta/2}. \qquad (3.48) $$

Hence, at least one of these two bounds holds whenever eq. (3.43) does.

<footnote n="3">³Where the 1/4 has been put just for convenience in order to alleviate some later expressions.</footnote>
We then turn to the second case in eq. (3.44), where obviously $\kappa \neq 0$, and we can hence rewrite eq. (3.38) as

$$ \gamma_{\pm} = 1 \pm \frac{\kappa}{|\kappa|} \frac{1}{\sqrt{1+\xi}}, \qquad (3.49) $$

with $\xi$ given by

$$ \xi = \frac{2\tilde{c}}{R^2\kappa^2} + \frac{\tilde{c}^2}{R^4\kappa^2}. \qquad (3.50) $$

The inequalities in eqs. (3.39) and (3.44) (and the assumption $\tilde{c} > 0$) imply in turn that<fnref n="4" />

$$ 0 < \xi \leq 2. \qquad (3.51) $$

Let us then first assume that $\kappa > 0$. In this case, we can obtain a useful lower bound on $\gamma_-$ noticing that, for $\xi$ in the range in eq. (3.51), one has

$$ 1 - \frac{1}{\sqrt{1+\xi}} \geq \left(1 - \frac{1}{\sqrt{3}}\right)^{\frac{\xi}{2}} > \frac{1}{8}\xi. \qquad (3.52) $$

Using this, we conclude that, for positive $\kappa$ and whenever eq. (3.44) holds, we have

$$ \gamma_- > \frac{\tilde{c}}{4R^2\kappa^2}. \qquad (3.53) $$

We also note, following the same logic, that $\gamma_+$ obeys exactly the same bound for negative $\kappa$ whenever eq. (3.44) holds. We can then, as previously, distinguish between the two subcases in eqs. (3.41) and (3.42) which implies that either $u^2 \geq |\kappa|R^2/2$ or $\tilde{v}^2 \geq |\kappa|R^2/2$. In the first subcase, and for positive $\kappa$, we can use the bounds in eqs. (3.29) and (3.53) to get that at large $R$

$$ \begin{aligned} \gamma_{-}f(u) &> \tilde{c}F_{0}2^{-\frac{\zeta}{2}}(\kappa R^{2})^{\frac{\zeta}{2}-1}\kappa^{-1} \\ &\geq \tilde{c}F_{0}2^{-\frac{\zeta}{2}}(R\sqrt{2\tilde{c}})^{\frac{\zeta}{2}-1}\kappa^{-1} \\ &\geq \tilde{c}F_{0}2^{-\frac{\zeta}{2}}(R\sqrt{2\tilde{c}})^{\frac{\zeta}{2}-1}, \end{aligned} $$

where to go from the first line to the second line we used once more eq. (3.44) and to go from the second line to the third line we used the fact that $1 \geq \kappa > 0$. Hence, using this last inequality as well as eq. (3.35), we get that, for positive $\kappa$ and whenever eq. (3.44) and $u^2 \geq \kappa R^2/2$ hold, we have

$$ \mathcal{U} > -2|g_0| + \tilde{c}F_0 2^{-\frac{\zeta}{2}} (R\sqrt{2\tilde{c}})^{\frac{\zeta}{2}-1}. \qquad (3.54) $$

Let us then turn to discuss the subcase where $\tilde{v}^2 > |\kappa|R^2/2$. In this subcase, still assuming that $\kappa > 0$ we notice that $\gamma_+ \geq 1$ and hence we get simply using eqs. (3.30) and (3.36) that

$$ \begin{aligned} \mathcal{U} &\geq -2|f_0| + 4G_0 2^{-\frac{\eta}{2}} (\kappa R^2)^{\frac{\eta}{2}} \\ &\geq -2|f_0| + 4G_0 2^{-\frac{\eta}{2}} (R\sqrt{2\tilde{c}})^{\frac{\eta}{2}}, \end{aligned} \qquad (3.55) $$

<footnote n="4">⁴This upper bound can be strengthened but is enough for the present discussion.</footnote>
where to go from the first line to the second line we used again eq. (3.44). Similar bounds can be obtained for negative $\kappa$, mutatis mutandis. Indeed, we find in this case, following the same logic, one has either

$$ U \geq -2|g_0| + 4F_0 2^{-\frac{\zeta}{2}} (R\sqrt{2\tilde{c}})^{\frac{\zeta}{2}} \quad (3.56) $$

or

$$ U > -2|f_0| + \tilde{c} G_0 2^{-\frac{\eta}{2}} (R\sqrt{2\tilde{c}})^{\frac{\eta}{2}-1}. \quad (3.57) $$

Hence, we conclude from the considerations of eqs. (3.47), (3.48) and (3.54) to (3.57) that the hypotheses (i) and (ii) above are enough to conclude that $\forall U_0 > 0, \exists R_0 > 0$ such that $R > R_0 \Rightarrow U > U_0$. This concludes our proof: in all cases, the “potential energy” $U$ grows without bound with $R$ implying that $|x|$ and $|y|$ stay bounded on shell.

To summarize, we have shown that the phase-space motion stays bounded at all times, irrespective of the initial conditions, for functions $f$ and $g$ which (i) are bounded below, and (ii) at large values of their arguments, obey respectively eq. (3.29) (where $F_0 > 0$ and $\zeta > 2$) and eq. (3.30) (where $G_0 > 0$ and $\eta > 2$). We stress that the behaviour of the functions $f$ and $g$, entering in condition (ii), only needs to be specified at large $|u|$ and $|v|$, implying, in particular, that the lower bounds $f_0$ and $g_0$ can be negative.

The above proof can also be extended to the case where, starting from the Liouville PP model, one ghostifies $x$ by applying a canonical transformation $\mathcal{C}_x^{\pm}$. The obtained models are in fact identical (up to a global sign in the Hamiltonian) to the previously considered PG models obtained by ghostifying $y$ from the PP model via a canonical transformation $\mathcal{C}_y^{\pm}$. As a result, the two kinds of models with a ghost are stable under the same conditions (i) and (ii). This is shown explicitly in appendix C.

## 3.2 A first application of the proof

As a first application of the above proof, we consider potentials $V_{LV}(x, y)$ generated by functions $f$ and $g$ which are polynomials in $u^2$ and $v^2$ respectively, i.e. of the form

$$ f(u) = \sum_{n=1}^{N_f} C_n (u^2)^n, \quad (3.58) $$

$$ g(v) = \sum_{n=1}^{N_g} \mathcal{D}_n (v^2)^n, \quad (3.59) $$

with $N_f \in \mathbb{N}$, $N_g \in \mathbb{N}$, $C_n \in \mathbb{R}$ and $\mathcal{D}_n \in \mathbb{R}$. More specifically, choosing $N_f = N_g = 2$, we get

$$ f(u) = \mathcal{C}_0 + \mathcal{C}_1 u^2 + \mathcal{C}_2 u^4, \quad (3.60) $$

$$ g(v) = \mathcal{D}_0 + \mathcal{D}_1 v^2 + \mathcal{D}_2 v^4, \quad (3.61) $$

and the potential $V_{LV}(x, y)$ is given by

$$ V_{LV}(x, y) = \frac{C_1 + \mathcal{D}_1}{2} + \frac{C_2 + \mathcal{D}_2}{2} W + \frac{C_0 - \mathcal{D}_0}{\tilde{W}} + (C_1 - \mathcal{D}_1) \frac{W}{2\tilde{W}} + \frac{C_2 - \mathcal{D}_2}{4} \left( \frac{W^2}{\tilde{W}} + \tilde{W} \right), \quad (3.62) $$
where $W$, $\tilde{W}$ and $w$ are defined in the previous subsection. Note that, in the above, the first term on the right-hand side is a mere constant and, as such, can be omitted. Moreover, the second term proportional to $W$ is a quadratic term for $x$ and $y$ (plus an irrelevant constant proportional to $\tilde{c}$). Choosing

$$C_2 > 0, \qquad (3.63)$$

$$D_2 > 0, \qquad (3.64)$$

$$\tilde{c} > 0, \qquad (3.65)$$

the functions $f$ and $g$ defined by eqs. (3.60) and (3.61) fulfill conditions (i) in eqs. (3.27) and (3.28) and (ii) in eqs. (3.29) and (3.30) of section 3.1 and as such yield models with bounded motion. The further choice $C_1 = D_1 = 1/2$, $C_2 = D_2 = 1/2$, $C_0 - D_0 = \lambda$, and $\tilde{c} = 1$ yields exactly the main model of reference [16] (whose Hamiltonian is given in equation (1) of this reference). The more general choice of a strictly negative $c < 0$, strictly positive and equal $C_2 = D_2 > 0$, and $C_1 \ge D_1$ (and, moreover, setting $D_1 = -C_1$ to remove an overall constant) gives the largest set of models introduced in [16] (i.e., in eq. (20) and (21) of this reference). The parameters $a, b, c$ and $d$ of [16] can be read off from the above as

$$a = -\frac{C_2 + D_2}{2},$$

$$b = C_0 - D_0,$$

$$c = -\frac{C_1 - D_1}{2},$$

$$d = \tilde{c},$$

where the left-hand side of the above refers to the notation of [16] (and we stress, in particular, that $c$ on the left-hand side of the third equality is introduced in [16] and is not the same as the $c$ of the present work). We see that the conditions eqs. (3.63) and (3.65) imply that $a < 0$ and $d > 0$, and the present work provides a stability proof that was not given in [16]. Note further that the condition on $c$ given in [16] (i.e. $c \le 0$) can, in fact, be relaxed and still yields stable motion.

### 3.3 A polynomial subclass of the Liouville model: towards more general conditions for stable motion

The Liouville model at the heart of this work, and defined in eqs. (3.2) and (3.5), possess an interesting subclass for which the potential $V_{LV}(x, y)$ is polynomial in $x$ and $y$. This subclass is generated by choosing

$$f(u) = \sum_{n=1}^{N} C_n (u^2)^n, \qquad (3.66)$$

$$g(v) = \sum_{n=1}^{N} C_n (v^2)^n, \qquad (3.67)$$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** We show the first two non-trivial examples of polynomial Liouville potentials, corresponding to $N = 3$ (left-hand panel) and $N = 4$ (right-hand panel) in eq. (3.68). With the notation of eq. (3.69): the left-hand panel, corresponds to $\omega_x^2 = 1$, $\omega_y^2 = -1$, $C_4 = 0$, and $\tilde{c} = 1$; the right-hand panel, corresponds to $\omega_x^2 = 1$, $\omega_y^2 = 1$, $C_4 = 1$, and $\tilde{c} = 1$. With regards to eq. (3.77), we also show $V_P(x) + V_G(y)$ as transparent surfaces. For clarity, we add respective contour plots of $(V_P + V_G) - V$ with contours at $(V_P + V_G) - V = (-8, -4, -2, 0, 2, 4, 8)$ and with respective shading from darker to lighter colour, where we also indicate regions in which $(V_P + V_G) - V \ge 0$.</figcaption></figure>

with $N \in \mathbb{N}$ (and still $C_n \in \mathbb{R}$), i.e., the same polynomial form for $f$ and $g$. In appendix D, we provide a proof that the resulting potential $V_{LV}(x, y)$, given by

$$V_{LV}^{(N)} = \sum_{n=1}^{N} \frac{C_n}{u^2 - v^2} \left[ (u^2)^n - (v^2)^n \right], \qquad (3.68)$$

is indeed polynomial in $x$ and $y$. To our knowledge, it was not noticed before that this particular choice leads to polynomials, and hence, for both the PP and PG cases, to integrable systems with polynomial interactions whose specific form can be found below and in the appendix D.

The explicit form of $V_{LV}(x, y)$, in the PG case, and for the simplest nontrivial case of interest here, is obtained from the expression in appendix D, where one trades the constants $C_2$ and $C_3$ for quadratic terms and chooses the constant $C_1$ to eliminate an irrelevant constant term. One obtains, ordering the interactions by powers of $x$ and $y$,

$$V_{LV}^{(4)}(x, y) = \frac{\omega_x^2}{2}x^2 - \frac{\omega_y^2}{2}y^2 + \frac{1}{\tilde{c}}\left(\frac{\omega_x^2}{2} - \frac{\omega_y^2}{2}\right)(x^2 - y^2)^2 + \tilde{c}C_4(x^4 - y^4) + C_4(x^2 - y^2)^3. \quad (3.69)$$

This rewriting captures the first two non-trivial cases, i.e., $N = 3$ and $N = 4$, for which we visualize examples of the potential in figure 3 and of the resulting phase-space motion in figure 4 (for $N = 3$) and figures 5 to 7 (for $N = 4$).

Except for the degenerate choice of $\omega_x = \omega_y$ (for which the $N = 3$ case reduces to the $N = 2$ case), the $N = 3$ case does not lead to fully stable motion, as conditions (i) in eqs. (3.27) and (3.28) and (ii) in eqs. (3.29) and (3.30) of section 3.1 are not fulfilled.<fnref n="5" /> Indeed, either one chooses $C_3$ positive and for large negative $v^2$, $g(v)$ diverges to $-\infty$, or one chooses $C_3$ negative and for large positive $u^2$, $f(u)$ diverges to $-\infty$. An example of two nearby sets of initial conditions — one stable and one unstable — is visualized in figure 4.

<footnote n="5">⁵Even though our proof in section 3.1 only provides sufficient not necessary conditions for the boundedness of motion, on physical grounds, we expect that integrals of motion unbounded from below and from above cannot bound the motion.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** We show trajectories (with initial conditions tuned to the boundary of an apparent island of stability) for $N = 3$, i.e., for $\omega_x^2 = 1$, $\omega_y^2 = -1$, $C_4 = 0$, and $\tilde{c} = 1$, in eq. (3.69). We superimpose a contour plot of the potential (cf. left-hand panel in figure 3) where increasingly light shades indicate increasingly large values of the potential. In the left-hand panel, we show a bounded trajectory. In the right-hand panel, we show an adjacent unstable trajectory. To be specific, and in case the reader wants to reproduce our results, we have chosen $x(t=0) = 1$, $x'(t=0) = 1$, $y(t=0) = A$, and $y'(t=0) = A$, for which the difference between bounded and unbounded motions is obtained by choosing respectively $A \approx 0.468$ in the left-hand panel and $A \approx 0.469$ in the right-hand panel. These values hence also locate the boundary between the two kinds of motions considered here.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** We show the trajectory (for random initial conditions drawn from a Gaussian distribution with standard deviation SD = 2) for $N = 4$, i.e., for $\omega_x^2 = 1$, $\omega_y^2 = 1$, $C_4 = 1$, and $\tilde{c} = 1$ in eq. (3.69). In the left-hand panel, we superimpose a contour plot of the potential (cf. right-hand panel in figure 3) where increasingly light shades indicate increasingly large (positive) values of the potential. In the right-hand panel, we show the respective evolution of the momenta.</figcaption></figure>

In contrast, the $N = 4$ case leads to bounded motion since, for $C_4 > 0$, conditions (i) in eqs. (3.27) and (3.28) and (ii) in eqs. (3.29) and (3.30) of section 3.1 are obeyed. We visualize three different examples: in figure 5, we show an example of the motion for equal frequencies, i.e., for $\omega_x^2 = 1$, $\omega_y^2 = 1$. In figure 6, we show an example of the motion for which the ghost has a tachyonic term, i.e., for $\omega_x^2 = 5$, but $\omega_y^2 = -5$. Finally, in figure 7, we show an example
<figure><img src="image_9.png" /><figcaption>**Figure 6.** As in figure 5 but for a case with a ghost frequency which is tachyonic around the origin, i.e., for $\omega_x^2 = 5$, $\omega_y^2 = -5$, $C_4 = 1$, and $\tilde{c} = 1$ in eq. (3.69).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** As in figure 5 but for a case with large frequency ratio between positive-energy and negative-energy mode, i.e., for $\omega_x^2 = 1$, $\omega_y^2 = 100$, $C_4 = 1$, and $\tilde{c} = 1$ in eq. (3.69).</figcaption></figure>

of the motion with a significant frequency hierarchy, i.e., for $\omega_x^2 = 1$ and $\omega_y^2 = 100$. In all cases, $C_4 = 1$. The “potential” in eq. (3.32) that stabilizes the motion is bounded from below and grows without bound at large $x$ and $y$ for $C_4 > 0$. It is given here by

$$ \mathcal{U}_{\text{LV}}^{(4)}(x, y) = \frac{\tilde{c}}{2}(\omega_x^2 x^2 + \omega_y^2 y^2) + \frac{1}{2}(x^4 - y^4)(\omega_x^2 - \omega_y^2) + C_4 \tilde{c} \left[ (x^2 + y^2)(x^2 - y^2)^2 + \tilde{c}(x^4 + y^4) \right], \quad (3.70) $$

where we rewrote eq. (D.12) using $\tilde{c} > 0$. Note that both the above $\mathcal{U}_{\text{LV}}^{(4)}$ and $V_{\text{LV}}^{(4)}(x, y)$ are symmetric with respect to separate reflections

$$ x \to -x, \quad \text{and} \quad y \to -y, \qquad (3.71) $$

as they depend only on $x$ and $y$ via their squares $x^2$ and $y^2$. The potential $V_{\text{LV}}^{(4)}(x, y)$ has
up to power six self-interactions and PG cross-couplings

$$ \frac{(\omega_y^2 - \omega_x^2)}{\tilde{c}} x^2 y^2, \quad \text{and} \quad 3C_4 x^2 y^2 (y^2 - x^2), \qquad (3.72) $$

between the usual DoF $x$ and the ghost $y$. These cross-coupling terms vanish when one of the coordinates is zero. However, elsewhere these cross-couplings cannot be neglected in comparison with the self-interactions. Without cross-couplings the positive energy degree of freedom $x$ has a self-interaction potential

$$ V_P(x) = \frac{\omega_x^2 x^2}{2} + \left[ \frac{(\omega_x^2 - \omega_y^2)}{2\tilde{c}} + C_4 \tilde{c} \right] x^4 + C_4 x^6, \qquad (3.73) $$

while the ghost self-interaction potential is

$$ V_G(y) = - \left[ \frac{\omega_y^2 y^2}{2} + \left[ \frac{(\omega_y^2 - \omega_x^2)}{2\tilde{c}} + C_4 \tilde{c} \right] y^4 + C_4 y^6 \right]. \qquad (3.74) $$

Without interactions one could flip the sign in front of the ghost's Hamiltonian, so that uncoupled ghost effectively experiences inverted self-interaction potential ($-V_G(y)$). Note that ($-V_G(y)$) can be obtained from ($V_P(x)$) by changing $x$ to $y$ everywhere, including indices.

The next non-trivial polynomial and stable PG theory is obtained for $N = 6$ (and positive $C_6$) and the corresponding potential is given in appendix D.

More generally, all polynomial theories with even $N$ (and positive $C_N$) are stable according to the above proof. This becomes apparent when looking at large $u^2 \gg |C_n|$ (as well as large $-v^2 \gg |C_n|$), for which

$$ \begin{aligned} f(u) &\rightarrow C_N (u^2)^N, \\ g(v) &\rightarrow (-1)^N C_N (-v^2)^N. \end{aligned} \qquad (3.75) $$

In writing the second expression, we arrange minus signs such as to account for $-v^2 > 0$.

For odd $N$, the proof conditions cannot be fulfilled since (for any choice of $C_N$) either $f(u)$ or $g(v)$ is unbounded below. While islands of stable initial conditions are possible (cf. left-hand panel in figure 4), we expect that there are always some initial conditions for which divergent behaviour remains (cf. right-hand panel in figure 4).

So, how do the proof conditions relate to the behaviour at large $x^2 \gg \tilde{c}$ and $x^2 \gg |C_n|$ (as well as large $y^2 \gg \tilde{c}$ and $y^2 \gg |C_n|$)? Taking both limits,

$$ V_{\text{LV}}^{(N)}(x, y) \rightarrow C_N (x^2 - y^2)^{N-1}. \qquad (3.76) $$

Together with the proof conditions, i.e., even $N$ and $C_N > 0$, this implies:

* The decoupled potential $V_P(x)$ (see eqs. (2.7) and (3.73)) is bounded below and unbounded above and thus serves as a stable potential for the decoupled motion of the $x$ mode.
* Similarly, $V_G(y)$ (see eqs. (2.7) and (3.74)) is bounded above and unbounded below and thus serves as a stable potential for the decoupled motion of the y mode.

* At sufficiently large $x^2$ and $y^2$, the interaction potential $V(x, y)$ is bounded by the decoupled potential, i.e.,

$$V(x, y) \leq V_P(x) + V_G(y) \quad \forall |x| \geq |y|. \qquad (3.77)$$

* As we will see, the special case of $|x| = |y|$, for which the full and the decoupled potential agree at sufficiently large $x^2$ and $y^2$, seems to require an additional criterion.

To elucidate the latter subtlety further, we take a look at another class of integrable polynomial potentials.

## 3.4 An unstable polynomial subclass: refining more general conditions for stable motion

To elucidate the possibility of stable motion in the Liouville class further, we can compare it to the motion in another class of integrable two-particle Hamiltonians.

As explained in section 2 real-valued PG Hamiltonians can sometimes also be obtained from complex-valued PP Hamiltonians. This turns out to be the case for another class of integrable two-particle Hamiltonians, closely related to the Liouville model, cf. appendix A. The details of this integrable model are discussed in appendix E and, following [49], we refer to it as class 3 from here on.

In fact, just as the integrable Liouville model, the integrable class 3 contains a polynomial subclass which is analogously obtained, cf. appendix E. Introducing frequencies, removing a constant shift (as in eq. (3.69)), and specifying to equal frequencies (i.e., $\omega_x = \omega_y = \omega$), the first few terms can be written as

$$V_{(\text{class 3})}^{(4)}(x, y) \stackrel{\omega_x=\omega_y}{=} \frac{\omega^2}{2}(x^2 - y^2) + C_4(x^2 - y^2)^3 - C_4 c(x+y)(x-y)(x \pm y)^2. \quad (3.78)$$

The general expression with unequal frequencies is given in appendix E. For $C_4 = 0$, the potential $V_{(\text{class 3})}^{(4)}$ is (depending on the other parameter choices) either bounded from below or bounded from above. More generally, just as for the Liouville polynomial subclass, each odd N case is thus unstable because one of the decoupled potentials is necessarily unstable.

For $C_4 > 0$, cf. also figure 8, the potential $V_{(\text{class 3})}^{(4)}$ fulfills all the conditions stated around eq. (3.77): both of the decoupled potentials are stable and the full potential is bounded by the decoupled potential. Nevertheless, we find that this case exhibits runaways with polynomial growth rate and is thus not stable. We thus conclude that the conditions stated around eq. (3.77) are not sufficient to guarantee the absence of runaway behaviour.

We first look at the case ($\omega = 0$), for which we observe that the runaway behaviour occurs either along the $x-y = 0$ or along the $x+y = 0$ direction. In fact, the potential in figure 8 exhibits sets of saddle points (marked with thick lines) extending out to infinity

<footnote n="6">⁶For $N = 4$, $|x| \geq |y|$ in (3.77) should be understood as $|x| \geq |y| + \frac{\omega_y^2 - \omega_x^2}{3\bar{c}C_4} \frac{1}{|x|+|y|}$, which agrees with $|x| \geq |y|$ either in the limit $|x| + |y| \to \infty$ or for $\omega_y^2 = \omega_x^2$.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** As in figure 3 but for the class 3 potential in eq. (3.78) with $\omega_x^2 = \omega_y^2 = 0$ and $C_4 = 1$. The thick lines indicate the directions along which unstable motion occurs. Note that the potential exhibits sets of saddle points along flat directions extending to infinity, cf. right-hand panel in figure 3.</figcaption></figure>

along flat directions of the potential. As a consequence of this, the local potential around $x \pm y = 0$ is no longer dominated by the respective decoupled potentials.

This can be seen from the Hessian matrix,<fnref n="7" /> which vanishes at $x = \pm y$, i.e.,

$$ \text{Hess}(V_{(\text{class } 3)}^{(4)}(x, y)) \stackrel{x=\pm y}{=} \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix}. \qquad (3.79) $$

In contrast, the Hessian matrix of the decoupled potential (at $x = \pm y$) is non-vanishing, i.e.,

$$ \text{Hess}(V_{(\text{class } 3)}^{(4)}(x, 0) + V_{(\text{class } 3)}^{(4)}(0, y)) \stackrel{x=\pm y}{=} 30 x^4 \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. \qquad (3.80) $$

Along the saddle surfaces, the local potential $V_{(\text{class } 3)}^{(4)}$ does, therefore, not resemble the decoupled potential. In fact, it does not resemble any decoupled potential. Put differently, there exist regions in which decoupled interactions no longer dominate the potential at large phase-space variables.

For $V_{(\text{class } 3)}^{(4)}$ with non-vanishing but equal frequencies (and still for $C_4 > 0$), the issue of runaways is lifted along one of the two directions. Without loss of generality, picking the ‘-’ sign in the last term in eq. (3.78), the $x = -y$ direction is stable but the $x = y$ direction remains unstable. Once more, this can be understood in terms of the Hessian matrix. For $x = -y$, the Hessian matrix now reads

$$ \text{Hess}(V_{(\text{class } 3)}^{(4)}(x, y)) \stackrel{x=-y}{=} (\omega^2 + 24 C_4 y^2) \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. \qquad (3.81) $$

$$ \text{Hess}(V_{(\text{class } 3)}^{(4)}(x, y)) \equiv \begin{pmatrix} \partial_x \partial_x V_{(\text{class } 3)}^{(4)} & \partial_x \partial_y V_{(\text{class } 3)}^{(4)} \\ \partial_y \partial_x V_{(\text{class } 3)}^{(4)} & \partial_y \partial_y V_{(\text{class } 3)}^{(4)} \end{pmatrix}. $$

<footnote n="7">⁷For the case at hand, the Hessian matrix is defined by</footnote>
This agrees with the Hessian matrix of a free theory, to be specific, the free theory with $V(x, y) = \frac{\omega^2}{2}(x^2 - y^2) + 2C_4(x^4 - y^4)$. In contrast, for $x = y$, we find

$$ \text{Hess}(V_{(\text{class 3})}^{(4)}(x, y)) \stackrel{x=y}{=} \omega^2 \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, \qquad (3.82) $$

for which the behaviour at large phase-space variables still has no free-theory equivalent. In agreement with the above, we find that, with equal frequencies, polynomial runaway behaviour no longer occurs along the $x = -y$ but still occurs along the $x = y$ direction.

Numerically, we observe that the runaways can be fully removed by the addition of decoupled interactions of the form $x^4 - y^4$. (Numerics are necessary since the addition of these terms breaks integrability.) With this addition, the Hessian matrix has the same structure as for the Liouville potential $V_{LV}^{(4)}(x, y)$. In both cases, it does still not agree with the Hessian matrix of the respective decoupled potential. However, one can identify a decoupled potential for which the Hessian matrices (at sufficiently large $x = \pm y$) agree. We conclude that it is sufficient (at least for the specific case at hand) that the Hessian matrix at large $|x| = |y|$ agrees with some decoupled potential, not necessarily with the decoupled potential of the respective model.

In passing, we mention that similar polynomial runaway behaviour also occurs in an integrable model previously obtained (from a supersymmetric construction) by Smilga et al. [21, 51]. Its Hamiltonian is given by

$$ H = \frac{1}{2}(p_x^2 + \omega^2 x^2) - \frac{1}{2}(p_y^2 + \omega^2 y^2) + \frac{\lambda}{4\omega}(x-y)(x+y)^3, \qquad (3.83) $$

where $\lambda$ and $\omega$ are constants. It follows as a degenerate limit of eq. (3.78), i.e., when picking the + sign in the last term of eq. (3.78) and taking $C_4 \to 0$ as well as $c \to \infty$ such that $C_4 c = -\lambda/4\omega$. (Alternatively, eq. (3.83) can also be obtained as a specific case<fnref n="8" /> of class 8 in appendix A.) The motion of this theory does not stay bounded (even though the ghost there is argued to be “benign” due to the slow (polynomial) growth-rate of the runaway solutions found numerically [21, 51]). In contrast to eq. (3.78), the potential of eq. (3.83) also violates the condition in eq. (3.77).

We conclude that all of the discussed cases can be understood if the following criterion is added to the list in section 3.3:

* For $|x| = |y|$ (or more generally, for all $(x, y)$ for which $V(x, y) = V_P(x) + V_G(y)$), the Hessian matrix at large phase-space variables approaches the Hessian matrix of some free theory.

Having understood these more subtle cases, we conclude: by translating the proof conditions to conditions on the potential $V(x, y)$, we understand that stable motion seems to be possible because stable decoupled potentials $V_P(x)$ and $V_G(y)$ dominate all interaction terms at sufficiently large $x^2$ and $y^2$.

<footnote n="8">⁸To be specific, the Hamiltonian in eq. (3.83) can be obtained from the integrable class 8 in eq. (A.15) by choosing $g = 0$ and $f'' = \frac{1}{2}\omega^2 + \frac{\lambda}{\omega}(x + iy)^2$. This results in the potential (for a PP theory) given by $V = (x^2 + y^2)(\frac{1}{2}\omega^2 + \frac{\lambda}{4\omega}(x + iy)^2)$. We then obtain the theory in eq. (3.83) by applying a canonical transformation $C_y^-$.</footnote>
In the following section 4, we will use this insight to formulate general conditions for stability. Before doing so, we comment on the relation to local Lyapunov stability.

## 3.5 Lagrange vs Lyapunov stability

In the previous subsections, we were concerned with global boundedness of motion which corresponds to the stability in the sense of Lagrange. However, one can be also interested in local stability in the sense of Lyapunov, see, e.g., [37]. In particular, this is crucial to specify the nature of equilibrium points which are fixed points of the system of equations of motion and can be considered as vacua. The two notions of stability are not equivalent. Indeed, a Lagrange stable system may not have locally stable solutions, while a system with Lyapunov stable equilibrium points may have other solutions with runaways invalidating stability in the sense of Lagrange. Here we will show that the Liouville class of systems in eq. (3.1), describing ghosts interacting with usual degrees of freedom ($\sigma = -1$), contains plenty of theories with Lyapunov stable equilibrium points. More specifically, we focus here first on systems in this class with functions $f(u)$ and $g(v)$ which are respectively analytic functions of $u^2$ and $\tilde{v}^2$ around the origin $(x, y) = (0, 0)$. These systems contain the polynomial case eqs. (3.66) and (3.67) as a subcase.

For such theories, the origin $(x, y) = (0, 0)$ is always an equilibrium point. Indeed, at the origin one has $u^2 = 0$ and $\tilde{v}^2 = \tilde{c}$, while

$$ \partial_x V_{LV} = 0, \quad \text{and} \quad \partial_y V_{LV} = 0, \qquad (3.84) $$

because all partial derivatives

$$ \frac{\partial u^2}{\partial x} = \left[ \frac{r^2 + \tilde{c}}{\tilde{W}} + 1 \right] x, \quad \frac{\partial \tilde{v}^2}{\partial x} = \left[ \frac{r^2 + \tilde{c}}{\tilde{W}} - 1 \right] x, \qquad (3.85) $$

$$ \frac{\partial u^2}{\partial y} = \left[ \frac{r^2 - \tilde{c}}{\tilde{W}} - 1 \right] y, \quad \frac{\partial \tilde{v}^2}{\partial y} = \left[ \frac{r^2 - \tilde{c}}{\tilde{W}} + 1 \right] y, \qquad (3.86) $$

are vanishing there. Now we can linearize the equations of motion<fnref n="9" />

$$ \ddot{x} + \frac{\partial V_{LV}}{\partial x} = 0, \quad \text{and} \quad \ddot{y} - \frac{\partial V_{LV}}{\partial y} = 0, \qquad (3.87) $$

around the origin and obtain corresponding frequencies

$$ \bar{\omega}_x^2 = \frac{2}{\tilde{c}^2} [\tilde{c}f'(0) - f(0) + g(\tilde{c})], \qquad (3.88) $$

$$ \bar{\omega}_y^2 = \frac{2}{\tilde{c}^2} [\tilde{c}g'(\tilde{c}) - g(\tilde{c}) + f(0)]. \qquad (3.89) $$

Note that $\bar{\omega}_x$ and $\bar{\omega}_y$ are just defined expanding around the origin and are hence conceptually different from unbarred frequencies defined previously, hence the different notation. Of course they could be the same depending on the potential. Naive stability by the linearized approximation requires that both $\bar{\omega}_x^2$ and $\bar{\omega}_y^2$ are positive. In particular, adding the two equations above, we see that stability requires that

$$ g'(\tilde{c}) + f'(0) > 0. \qquad (3.90) $$

<footnote n="9">⁹Note that the sign difference in the y-equation of motion appears due to the negative mass of the ghost.</footnote>
Let us now confirm that the conditions $\bar{\omega}_x^2 > 0$ and $\bar{\omega}_y^2 > 0$ are sufficient for the Lyapunov stability of the equilibrium point in the origin. Consider the first integral $J_{LV}$ given by eq. (3.31). Expanding $J_{LV}$ around the origin and around vanishing momenta one obtains

$$J_{LV} = g(\tilde{c}) + f(0) + \frac{\tilde{c}}{2}(p_x^2 + p_y^2) + \frac{\bar{\omega}_y^2 \tilde{c}}{2} y^2 + \frac{\bar{\omega}_x^2 \tilde{c}}{2} x^2 + \dots, \quad (3.91)$$

where ellipsis stands for higher order terms. Clearly, for $\bar{\omega}_x^2 > 0$ and $\bar{\omega}_y^2 > 0$ the first integral has a minimum at the equilibrium point in the origin. Thus, $J_{LV} - (g(\tilde{c}) + f(0))$ is a Lyapunov function. Indeed, it is conserved, vanishing at the equilibrium solution and is positive definite in an open neighbourhood around it (see [37]). Thus, the origin is Lyapunov stable. Here, the conserved quantity $J_{LV}$ plays the analog of usual energy for establishing stability. Crucially, the “kinetic energies” are all positive definite in $J_{LV}$. Moreover, $U(x, y)$ — the “potential part” of the “energy” — has an isolated minimum at the origin. Thus, exchanging the usual potential energy with $U(x, y)$, one meets the conditions of the classical Lagrange theorem stating that the minimum of the potential energy corresponds to a stable equilibrium configuration. In contrast, it is interesting to note that this stable equilibrium point is a saddle point for the potential $V_{LV}(x, y)$ as

$$V_{LV} = \frac{f(0) - g(\tilde{c})}{\tilde{c}} + \frac{\bar{\omega}_x^2 x^2}{2} - \frac{\bar{\omega}_y^2 y^2}{2} + \dots \quad (3.92)$$

It also appears above that $\bar{\omega}_x^2 > 0$ and $\bar{\omega}_y^2 > 0$ do guarantee Lyapunov stability for the usual PP systems due to the Lagrange theorem. The Lyapunov stability of the origin, implies, by continuity, the same stability for solutions in a small neighbourhood. It is important to repeat that the conditions on $f(u)$ and $g(v)$ required for the Lyapunov stability of the origin (positivity of (3.88) and (3.89)) can be easily violated in systems which are globally stable according to Lagrange. And vice versa: a Lagrange unstable system may have a Lyapunov stable origin. This can be well illustrated by the case of $V_{LV}^{(3)}$ PG system given by eq. (3.68), see also eq. (D.4). As we discussed before, this system is not stable from the point of view of Lagrange, however the origin is Lyapunov stable provided $\omega_x^2$ and $\omega_y^2$ chosen in eq. (D.6) are positive. On the other hand, as we discussed before, $V_{LV}^{(4)}$ from eq. (3.69) is Lagrange stable for $C_4 > 0$ and $\tilde{c} > 0$, but can have a Lyapunov unstable origin if one of $\omega_x^2, \omega_y^2$ chosen in eq. (D.6) is negative. To finish, let us here mention that the Lyapunov stability of the origin equilibrium point for a distinct class of interacting ghosty systems has previously been demonstrated in [28] using different methods. Moreover, in [16] the Lyapunov stability of the origin has been proven for a particular PG system from the same Liouville class in eq. (3.1) using methods similar to those employed in the current paper.

**All vacua and their Lyapunov stability for the polynomial subclass of the Liouville model.** The PG interacting system with polynomial potential in eqs. (3.69) and (D.7) belongs to the class discussed above and deserves to be analyzed in more detail. In figure 5, figure 6 and figure 7 we have already presented numerical solutions of this system for various parameters. In this subsection, we find all equilibrium points, or fixed points of the evolution $(x_e, y_e)$, i.e.,

$$\partial_x V_{LV}^{(4)}(x_e, y_e) = \partial_y V_{LV}^{(4)}(x_e, y_e) = 0. \quad (3.93)$$
We then determine whether these vacua are Lyapunov stable. Assuming non vanishing $C_4$ and $\tilde{c}$, the above two equations are equivalent to the vanishing of

$$x_e \left[ x_e^2 + \frac{(x_e^2 - y_e^2)}{2\tilde{c}} \left[ \frac{(\omega_x^2 - \omega_y^2)}{C_4\tilde{c}} + 3(x_e^2 - y_e^2) \right] + \frac{\omega_x^2}{4C_4\tilde{c}} \right],$$ $$y_e \left[ y_e^2 + \frac{(x_e^2 - y_e^2)}{2\tilde{c}} \left[ \frac{(\omega_x^2 - \omega_y^2)}{C_4\tilde{c}} + 3(x_e^2 - y_e^2) \right] + \frac{\omega_y^2}{4C_4\tilde{c}} \right],$$

whose solutions will be discussed below. We then shall demonstrate how PG interactions which are supposed to destabilize the system, instead introduce new stable vacua. Interestingly, the appearance of these new vacua leads to spontaneous breaking of the reflection symmetry in eq. (3.71). To establish Lyapunov stability, we will use the result of appendix B: the extrema of $V_{LV}^{(4)}$ are also extrema of $U_{LV}^{(4)}$, see eq. (B.4). Then, similarly to the discussion above on the stability of the origin, the minimum of $U_{LV}^{(4)}$ implies the Lyapunov stability of the vacuum, as $J_{LV}$ from eq. (3.31) can be taken again as a Lyapunov function.

The first equilibrium point is the origin $(x_e, y_e) = (0, 0)$, which was discussed at the beginning of this section, and the Lyapunov stability just requires that $\omega_x^2 > 0$ and $\omega_y^2 > 0$ which we shall also assume to hold below. Note that this vacuum respects the reflection symmetry in eq. (3.71) for both coordinates and is such that (see eq. (3.70))

$$U_{LV}^{(4)}(x_e, y_e) = 0. \qquad (3.94)$$

One can find a second type of equilibrium point where only one of the coordinates vanishes but not the other one. Indeed, $x_e = 0$, or $y_e = 0$, is obviously a solution to eq. (3.93) (see the two above equations). When either $x_e = 0$, or $y_e = 0$, the stability along the direction where the equilibrium coordinate does not vanish can be analyzed by looking at the potential in eq. (3.74) or eq. (3.73). It is seen there that, for sufficiently small $C_4$, either the coefficient in front of $x^4$ in $V_P(x)$ or in front of $y^4$ in $(-V_G(y))$ can become negative, indicating that the origin is not the only minimum of $V_P(x)$ or of $(-V_G(y))$. In that case, the vacuum spontaneously breaks the reflection symmetry in eq. (3.71) with respect to this coordinate. Let us assume for concreteness that $\omega_x < \omega_y$, so that we take $y_e = 0$ as then implied by the positivity of the coefficient in front of $y^4$ in eq. (3.74). The opposite case is obtained from all formulas with the replacement $x \leftrightarrow y$. Looking for non-vanishing solutions of $\partial_x V_P(x) = 0$, yields the following quartic equation

$$\omega_x^2 + 2 \left[ \frac{\omega_x^2 - \omega_y^2}{\tilde{c}} + 2C_4 \tilde{c} \right] x^2 + 6C_4 x^4 = 0, \qquad (3.95)$$

which has real non-vanishing solutions $x^2$ for

$$\left[ \frac{\omega_x^2 - \omega_y^2}{\tilde{c}} + 2C_4 \tilde{c} \right]^2 > 6C_4 \omega_x^2. \qquad (3.96)$$

The last condition is easy to fulfill for a sufficiently small $C_4$, i.e., for (see also footnote 11)

$$4\tilde{c}^2 C_4 < \omega_x^2 + 2\omega_y^2 - \omega_x \sqrt{3(4\omega_y^2 - \omega_x^2)}. \qquad (3.97)$$
The solution of eq. (3.95) is

$$ \frac{x_{b(\pm)}^2}{\tilde{c}} = \frac{\Omega_y^2 - \Omega_x^2 - 1}{3} \pm \sqrt{\left(\frac{\Omega_y^2 - \Omega_x^2 - 1}{3}\right)^2 - \frac{\Omega_x^2}{3}}, \quad (3.98) $$

where we denoted<fnref n="10" />

$$ \Omega_x^2 = \frac{\omega_x^2}{2C_4\tilde{c}^2}, \quad \text{and} \quad \Omega_y^2 = \frac{\omega_y^2}{2C_4\tilde{c}^2}. \quad (3.99) $$

To obtain a real $x_{b(\pm)}$, one has to require that the first term in the right-hand side of eq. (3.98) is positive<fnref n="11" /> which translates to

$$ 2C_4\tilde{c}^2 < \omega_y^2 - \omega_x^2. \quad (3.100) $$

Note that this inequality in eq. (3.100) is weaker than eq. (3.97), so that, when the latter holds, both extrema $x_{b(\pm)}^2$ are real and in total there are four solutions. Clearly, $x_{b(-)}^2$ corresponds to a local maximum of $V_P(x)$ while $x_{b(+)}^2$ corresponds to a local minimum and hence a stable direction along the $x \neq 0$ direction.

Now we can investigate when such points $(\pm x_{b(+)}, 0)$ are minima of $\mathcal{U}_{LV}^{(4)}$ given by eq. (3.70). First one notices, that at these points $\partial_x \partial_y \mathcal{U}_{LV}^{(4)}(\pm x_{b(+)}, 0) = 0$, while the other second derivatives are

$$ \partial_x^2 \mathcal{U}_{LV}^{(4)} = \frac{8}{3} C_4 \tilde{c}^3 - \frac{4}{3} \tilde{c} (\omega_x^2 + 2\omega_y^2) + \frac{2(\omega_x^2 - \omega_y^2)^2}{3C_4 \tilde{c}} - \frac{2}{3} \mathcal{R} \left[ 2\tilde{c} + \frac{\omega_x^2 - \omega_y^2}{C_4 \tilde{c}} \right], $$

and

$$ \partial_y^2 \mathcal{U}_{LV}^{(4)} = \frac{1}{2} \tilde{c} \left[ 2\omega_y^2 - \frac{(\mathcal{R} - 2C_4\tilde{c}^2 - \omega_x^2 + \omega_y^2)^2}{9C_4\tilde{c}^2} \right], $$

where we denoted

$$ \mathcal{R} = \sqrt{4C_4^2 \tilde{c}^4 + (\omega_x^2 - \omega_y^2)^2 - 2C_4(\omega_x^2 + 2\omega_y^2)\tilde{c}^2}. $$

One can deduce that both second derivatives are positive provided

$$ \frac{(\omega_x^2 - \omega_y^2)^2}{2\omega_y^2} < 4C_4\tilde{c}^2 < \omega_x^2 + 2\omega_y^2 - \omega_x\sqrt{3(4\omega_y^2 - \omega_x^2)}. \quad (3.101) $$

Note that for the assumed $\omega_y > \omega_x$ the lower bound above is always smaller than the upper bound. Thus there is always a space for $C_4$ between them.<fnref n="12" /> Moreover, the upper bound is

<footnote n="10">¹⁰This form reveals the essence of rescalings in eq. (3.8).</footnote>

<footnote n="11">¹¹Note that this inequality forbids $4\tilde{c}^2C_4 > \omega_x^2 + 2\omega_y^2 + \omega_x\sqrt{3(4\omega_y^2 - \omega_x^2)}$ which would be another way to satisfy eq. (3.96).</footnote>

<footnote n="12">¹²Note that one can reinterpret this inequality as $(\Omega_x^2 - \Omega_y^2)^2/(2\Omega_y^2) < 2 < \Omega_x^2 + 2\Omega_y^2 - \Omega_x\sqrt{3(4\Omega_y^2 - \Omega_x^2)}$ for rescaled frequencies in eq. (3.99), see figure 9.</footnote>
<table><caption>**Table 1.** We list all possible minima of $\mathcal{U}_{\text{LV}}^{(4)}(x, y)$, from eq. (3.70) where we assumed that $\tilde{c} > 0$, $C_4 > 0$ and $\Omega_{x/y}^2 > 0$. In the last column, we list conditions for the existence of such minima. These minima are Lyapunov stable equilibrium configurations or stable vacua of the theory. The origin respects the reflection symmetry in eq. (3.71). The “decoupled-$x$” minima are given by eq. (3.98) and break the reflection symmetry $x \to -x$. The “decoupled-$y$” minima are given by eq. (3.98), with exchange $x \leftrightarrow y$, and break the reflection symmetry $y \to -y$. The “coupled” minima are given by eqs. (3.102) and (3.103) and break both of the reflection symmetries in eq. (3.71).</caption><thead><tr><th></th><th>x<sup>2</sup></th><th>y<sup>2</sup></th><th>conditions</th></tr></thead><tbody><tr><td>origin</td><td>0</td><td>0</td><td>&Omega;<sub>x/y</sub><sup>2</sup> &gt; 0,</td></tr><tr><td>decoupled-x</td><td><sup>c&#773;</sup>&frasl;<sub>3</sub> [(&Omega;<sub>y</sub><sup>2</sup> &minus; &Omega;<sub>x</sub><sup>2</sup> &minus; 1) + &radic;(&Omega;<sub>y</sub><sup>2</sup> &minus; &Omega;<sub>x</sub><sup>2</sup> &minus; 1)<sup>2</sup> &minus; 3&Omega;<sub>x</sub><sup>2</sup>]</td><td>0</td><td>(3.101)</td></tr><tr><td>decoupled-y</td><td>0</td><td><sup>c&#773;</sup>&frasl;<sub>3</sub> [(&Omega;<sub>x</sub><sup>2</sup> &minus; &Omega;<sub>y</sub><sup>2</sup> &minus; 1) + &radic;(&Omega;<sub>x</sub><sup>2</sup> &minus; &Omega;<sub>y</sub><sup>2</sup> &minus; 1)<sup>2</sup> &minus; 3&Omega;<sub>y</sub><sup>2</sup>]</td><td>(3.101) x &harr; y</td></tr><tr><td>coupled</td><td><sup>c&#773;</sup>&frasl;<sub>8</sub> [(&Omega;<sub>x</sub><sup>2</sup> &minus; &Omega;<sub>y</sub><sup>2</sup>)<sup>2</sup> &minus; 4&Omega;<sub>x</sub><sup>2</sup>]</td><td><sup>c&#773;</sup>&frasl;<sub>8</sub> [(&Omega;<sub>x</sub><sup>2</sup> &minus; &Omega;<sub>y</sub><sup>2</sup>)<sup>2</sup> &minus; 4&Omega;<sub>y</sub><sup>2</sup>]</td><td>(3.104)</td></tr></tbody></table>

just what is required for the very existence of the equilibrium point, see eq. (3.97). Hence, we have just shown that, for $C_4$ satisfying the bound above, the stable points $(\pm x_{b(+)}, 0)$ correspond to local minima of $\mathcal{U}_{\text{LV}}^{(4)}$ (see figure 11) and these vacua are Lyapunov stable, as they actually also are turning off interactions with the ghost $y$. For the opposite hierarchy between frequencies, analogous equilibrium points are $(0, \pm y_{b(+)})$ and the whole analysis and formulas are applicable after exchange of $x \leftrightarrow y$. Interestingly, for $C_4$ smaller than the lower bound in eq. (3.101), the equilibrium points $(\pm x_{b(+)}, 0)$ are not local minima of $\mathcal{U}_{\text{LV}}^{(4)}$ so that one cannot warrant stability. This is a strong hint that these equilibrium configurations then become unstable. On the other hand, one can check that $(\pm x_{b(-)}, 0)$ never become minima of $\mathcal{U}_{\text{LV}}^{(4)}$ showing that these equilibrium configurations remain always unstable. Of course the same logic is applicable after exchange of $x \leftrightarrow y$.

Finally, the last category of equilibrium points appears for such small coupling constant $C_4$ that the PG-interactions in eq. (3.72) introduce new nontrivial equilibrium points. These vacua spontaneously break both reflections in eq. (3.71) and are solutions of the system in eq. (3.93)

$$x^2 + \frac{(\omega_x^2 - \omega_y^2)(x^2 - y^2)}{2C_4\tilde{c}^2} + 3\frac{(x^2 - y^2)^2}{2\tilde{c}} + \frac{\omega_x^2}{4C_4\tilde{c}} = 0,$$

$$y^2 + \frac{(\omega_x^2 - \omega_y^2)(x^2 - y^2)}{2C_4\tilde{c}^2} + 3\frac{(x^2 - y^2)^2}{2\tilde{c}} + \frac{\omega_y^2}{4C_4\tilde{c}} = 0.$$

By subtracting these equations from each other one obtains

$$4C_4\tilde{c}(y^2 - x^2) = \omega_x^2 - \omega_y^2,$$

so that now one can substitute $y^2 - x^2$ back into both equations to obtain, using the same
<figure><img src="image_12.png" /><figcaption>**Figure 9.** In the left-hand panel, using the rescaled frequencies in eq. (3.99), we visualize the conditions for the existence of different stable vacua listed in table 1. In the white region, the origin is the only Lyapunov-stable vacuum. In the yellow stripe above (below), the stable origin is supplemented by the two Lyapunov-stable “decoupled-x” (“decoupled-y”) vacua. Further up and/or to the right, in the orange region, the stable origin is supplemented by the four Lyapunov-stable coupled vacua. In the other three panels, we specify to $C_4 = 1$ as well as $\tilde{c} = 1$ and show the three indicated examples for the respective possible cases. We super-impose the latter on contour plots of $\log(|\mathcal{U}_{\text{LV}}^{(4)}(x, y)|)$. (The contours are spaced by half-integers with the thick-dashed contour indicating $\log(|\mathcal{U}_{\text{LV}}^{(4)}(x, y)| = 3$.) We omit further configurations which arise from exchanging $(x \leftrightarrow y)$. From left to right, we show $(\Omega_x^2, \Omega_y^2) = (0.5, 7.5)$ (i.e., the point labeled “I” in the orange region of the left-hand panel), $(\Omega_x^2, \Omega_y^2) = (0.5, 3.75)$ (i.e., the point labeled “II” in the yellow region of the left-hand panel) and $(\Omega_x^2, \Omega_y^2) = (0.5, 0.5)$ (i.e., the point labeled “III” in the white region of the left-hand panel). Larger (smaller) dots indicate minima (other real-valued extrema) of $\mathcal{U}_{\text{LV}}^{(4)}(x, y)$ and the colours white, yellow, and orange indicate origin, decoupled extrema, and coupled extrema, respectively. Even though, for illustrative purposes, we have not chosen the parameters to preserve the ratio of $\Omega$'s but rather to fit the minima at the same region of $(x, y)$ space, it is clear that, for fixed values of $\tilde{c}$ and $\omega_y$, $\omega_x$, the vacuum configurations will change from left to right with increasing strength of $C_4$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** For values of $C_4$ satisfying eq. (3.104), we show a logarithmic plot of the “uplifted” effective “potential” $\mathcal{U}_{\text{LV}}^{(4)}(x, y)$, i.e., the plot of $\log(1 + \mathcal{U}_{\text{LV}}^{(4)}(x, y) - \mathcal{U}_{\text{LV}}^{(4)}(x_I, y_I))$, such that we take the logarithm of a strictly positive quantity. The red points indicate the four global vacua $(x_I, y_I)$ given by eqs. (3.102) and (3.103) with all combinations of signs for the square roots. The black point indicates the vacuum at the origin which, however, is only a local minimum of $\mathcal{U}_{\text{LV}}^{(4)}$. The parameters are chosen as $\tilde{c} = 1$, $\omega_x = 1$, and $\omega_y = 1.7$, while $C_4$ is 1.5 times smaller than the lower bound in eq. (3.104) so that $C_4 \approx 0.103$. This case has different numerical values of constants, but its structure of vacua corresponds to case I in figure 9.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 11. On the left, we show the effective “potential” $U_{LV}^{(4)}(x, y)$ around vacua, for $C_4$ located in the interval (3.101), while on the right, we plot, in the logarithmic scales, the same “uplifted” effective “potential”: $\log(1 + U_{LV}^{(4)}(x, y) - U_{LV}^{(4)}(x_{b(+)}, 0))$, on larger scales. The red points indicate the two global vacua $(\pm x_{b(+)}, 0)$ from eq. (3.98). The black point indicates the trivial vacuum at the origin, which however is only a local minimum of $U_{LV}^{(4)}$. The parameters are set to $\tilde{c} = 1$, $\omega_x = 1$, and $\omega_y = 1.7$, as in figure 10, while $C_4$ is higher than that of figure 10 and is in the middle of the interval (3.101): $C_4 \approx 0.221$. This case has different numerical values of constants, but its structure of vacua corresponds to case II in figure 9.</figcaption></figure>

definitions as in eq. (3.99),

$$ \frac{x_I^2}{\tilde{c}} = \frac{(\Omega_x^2 - \Omega_y^2)^2 - 4\Omega_x^2}{8}, \qquad (3.102) $$

and similarly

$$ \frac{y_I^2}{\tilde{c}} = \frac{(\Omega_x^2 - \Omega_y^2)^2 - 4\Omega_y^2}{8}. \qquad (3.103) $$

The right-hand side of both equations are strictly positive provided that<fnref n="13" /> (as easily seen using eq. (3.99))

$$ C_4 < \min \left( \frac{(\omega_x^2 - \omega_y^2)^2}{8\omega_x^2\tilde{c}^2}, \frac{(\omega_x^2 - \omega_y^2)^2}{8\omega_y^2\tilde{c}^2} \right), \qquad (3.104) $$

which guarantees that these nontrivial equilibrium points exist. Further, this inequality also implies that $\partial_x^2 U_{LV}^{(4)} > 0$ along with the determinant of the Hessian $H_{css}(U_{LV}^{(4)}) > 0$, see eq. (B.13). Thus, all the four equilibrium points<fnref n="14" /> $(x_I, y_I)$, $(-x_I, y_I)$, $(x_I, -y_I)$, $(-x_I, -y_I)$

<footnote n="13">¹³These inequalities trivially transform to inequalities on rescaled frequencies, e.g., $(\Omega_x^2 - \Omega_y^2)^2 > 4\Omega_y^2$, see figure 9.</footnote>

<footnote n="14">¹⁴Here we assumed that the square root is taken as a principal square root returning non-negative numbers.</footnote>
are Lyapunov stable, as the “potential” $\mathcal{U}_{LV}^{(4)}$ has local minima there. One can check that

$$12\mathcal{C}_4^3\tilde{c}^4 \mathcal{U}_{LV}^{(4)}(x_I, y_I) = -\left[32\mathcal{C}_4^2(\omega_x^4 + \omega_y^4)\tilde{c}^4 + -8\mathcal{C}_4(\omega_x^2 - \omega_y^2)^2(\omega_x^2 + \omega_y^2)\tilde{c}^2 + (\omega_x^2 - \omega_y^2)^4\right] < 0.$$

Thus, comparing the above to eq. (3.94), we see that the equilibrium points above are in a sense true vacua as $\mathcal{U}_{LV}^{(4)}(x_I, y_I) < \mathcal{U}_{LV}^{(4)}(x_e, y_e)$, see also figure 10.

We would like to stress that we found all possible equilibrium points, i.e. solutions to eq. (3.93). Moreover, comparing, eq. (3.101) with eq. (3.104) one infers that the existence of stable equilibrium points $(\pm x_I, \pm y_I)$ and $(\pm x_{b(+)}, 0)$ (or $(0, \pm y_{b(+)})$) is mutually exclusive. Hence, PG-interactions in eq. (3.72) can actually introduce new Lyapunov-stable equilibrium points which would be impossible if the ghost $y$ had not interacted with the positive-energy DoF $x$. Thus, at least for discrete dynamical degrees of freedom, interactions with ghost can create new stable vacua. Finally, we would like to mention that it is the strength of the coupling constant $\mathcal{C}_4$ which controls the symmetry breaking pattern for fixed $\tilde{c}$ and $\omega_{x/y}$. The larger the constant, the more symmetry is there, see figure 9.

## 4 Numerical investigation of nonintegrable models

We now relax the technical assumption of integrability. In sections 3.3 and 3.4, we have rephrased the conditions for stability of the proof of section 3.1 (which are stated in terms of free functions, specific to the integrable model) to conditions in terms of the shape of the potential $V(x, y)$.

This motivates us to investigate whether similar conditions lead more generally to stable motion and, in particular, when there is no integrability.

We find numerically that: $(\aleph)$, a sufficient condition for the stable motion of a Hamiltonian point-particle system with coupled positive- and negative-sign kinetic terms for the respective coordinates $x$ and $y$ is that its full potential $V(x, y)$, at large coordinates $x$ and $y$, is sufficiently dominated everywhere by a potential of the form $\tilde{V}_P(x) + \tilde{V}_G(y)$ where both $\tilde{V}_P$ and $-\tilde{V}_G$ are yielding by themselves stable motions along $x$ and $y$. I.e. we need the interactions to be sufficiently suppressed with respect to a stable “decoupled” part at large values of the dynamical variables.

Note that these conditions are reminiscent of the Kolmogorov-Arnold-Moser (KAM) theorem (for a pedagogical discussion see [39, appendix 8]), as the theory with the Hamiltonian just given by $\tilde{V}_P(x) + \tilde{V}_G(y)$, which dominates at large distance, would be integrable (as it has two separately conserved quantities: the respective energies of $x$ and $y$), our statement is just that, it is enough for the interactions between $x$ and $y$ to be sufficiently subdominant at large $x$ and $y$ to ensure stability. The conditions stated in sections 3.3 and 3.4 are an explicit expression of the above statement for the specific integrable polynomial subclasses at hand. Beyond integrability, we still find that, whenever the above statement $(\aleph)$ holds, all (at least, all numerically investigated) initial conditions avoid runaways (at least, up to the finite numerical evolution time).
Naturally, the lack of integrability also means that the analytical proof strategy in section 3 will no longer apply. Hence, we need to resort to numerical evidence. The latter comes with the difficulty that — in the absence of rigorous analytical proof — runaways may only be encountered at very late evolution times and/or only for a small subset of the initial conditions. In section 4.1, we introduce a specific numerical setup to alleviate these issues, at least to some degree. In section 4.2, we demonstrate that the numerics are able to confidently distinguish between stable and unstable potentials in the presented subclass of integrable and polynomial Liouville models. With this benchmark at hand, we apply the same numerical analysis to several nonintegrable models in section 4.3, i.e., to two classes of non-polynomial and to one class of polynomial potentials. This provides evidence that the above statement (ical) holds true, even in the absence of integrability.

We emphasize that all of our results suggest that the statement (ical) only needs to hold outside of a finite region in phase space. As a special case, it thus encompasses all interaction potentials which are localised in a finite phase-space region and supplemented by standard quadratic terms $V_P = \frac{1}{2}\omega_x^2 x^2$ and $V_G = -\frac{1}{2}\omega_y^2 y^2$ with $\omega_x, \omega_y \in \mathbb{R}_+$. Two explicit examples of this are presented in section 4.3. Note that this contrasts with the usual conditions used in the KAM theorem which uses perturbations of an integrable system which are everywhere small with respect to the integrable part, and not only at large values of the dynamical variables.

We also highlight that, if the above holds also for nonintegrable systems, it would imply that the addition of sufficiently strong self-interactions can stabilize any PG system. Indeed, all of our numerics support this conclusion: an explicit example is discussed in section 4.3.

At the same time, we caution that there is a — potentially crucial — distinction between integrable and nonintegrable models. Integrability ensures that whenever the motion returns to the initial position $(x_0, y_0)$, also the momenta must return to their initial values. For nonintegrable systems, this need not be the case. Nonintegrable systems may therefore be distinct in that they can potentially extract a pair of positive and negative kinetic energy from the finite region in phase space in which the above conditions are violated. This may lead to a slow growth of the explored phase-space region. While we do not observe any evidence for such slow growth in our simulations, we cannot exclude it. In order to avoid such a growth, one could for example require that the Hamiltonian in the finite region be sufficiently close to integrable so that when the system returns to the initial position $(x_0, y_0)$, the momenta return to values sufficiently close to their initial values. Intriguingly, this does not a priori forbid arbitrarily strong interactions between a ghost and a positive energy degree of freedom in the finite region. Whether it is possible to make this kind of statement more quantitative and precise is an interesting question which could be investigated in the future.

## 4.1 Numerical setup

For our numerical experiments, we draw random initial conditions $x(t_i)$, $y(t_i)$, $x'(t_i)$, $y'(t_i)$ at $t_i = 0$ from a normal distribution with standard deviation SD. We then evolve the respective initial value problem with the Runge-Kutta method$^{15}$ up to either $t_f = t_{max}$ or until the

<footnote>$^{15}$We use the native NDSolve routine in Mathematica, specifying to an embedded Runge-Kutta 5(4) pair [52] following [53]. Further, we specify MaxSteps $\to \infty$ such that the integration is terminated either at $t = t_{max}$ or if NDSolve detects stiffness.</footnote>
numerical routine detects a runaway (see conditions below), in which case $t_f < t_{max}$. We make sure that $t_{max}$ is sufficiently large such as to differentiate between cases with and without runaways. We then repeat the experiment $N_{evol}$ times and determine the mean evolution time $t_{mean} = \frac{1}{N_{evol}} \sum_{j=1}^{N_{evol}} t_j^{(j)}$. Assuming that the criteria below confidently detect a runaway, a value of $t_{mean} = t_{max}$ signals the absence of any runaway solutions in the numerical experiment.

In the above numerical experiment, we distinguish between two criteria to detect runaways:

(A) We use stiffness detection of our numerical evolution scheme as a proxy for catastrophic (i.e., faster than exponential) runaways. Stiff numerical evolution typically occurs whenever very different scales enter the evolution. In practice, given that we evolve with a Runge-Kutta method, stiffness is detected by the dominant eigenvalue of the Jacobian, see [53, 54] and references therein.

(B) We find that polynomial growth is harder to detect since the eigenvalues of the Jacobian can remain bounded throughout the runaway. In particular, this also includes asymptotically flat potentials in which trajectories can escape to infinity once they escape the potential in a confined phase-space region. To detect such cases, we perform the same type of experiment but we add an additional abort criterion that detects if $x(t) > C$ or $y(t) > C$ outgrow a specified bound $C$. This criterion thereby detects whenever the evolution escapes a specified phase-space region. When we pick $C \gg$ SD, we expect that only rare cases of stable motion are falsely detected as unstable.

We will highlight below in which of the cases it makes a difference whether criterion (B) is applied in addition to criterion (A).

## 4.2 Benchmark: the Liouville model revisited

Here, we use the previously discussed integrable Liouville model, cf. section 3, to demonstrate that the numerical experiment meaningfully distinguishes between the presence and the absence of runaways.

We first focus on the polynomial subclass, cf. eqs. (3.68) and (3.75), for which the resulting potential $V(x, y)$ reduces to a polynomial of degree $N$. As discussed in section 3.3, the proof conditions imply stability at even degree $N$ (assuming also that $C_N > 0$) and, moreover, align with conditions stated in terms of the potential.

For each polynomial degree $2 \le N \le 9$, we perform our numerical experiment $N_{evol} = 10^2$ times. We choose a maximal evolution time $t_{max} = 10^2$ and draw random initial conditions from a normal distribution with standard deviation $SD = 1$.

The result is shown in figure 12 and confirms stability (instability) at even (odd) degree $N$. It turns out to be irrelevant whether or not we apply criterion (B). We take this as an indication that the runaways occur with faster than exponential growth.

Going beyond the polynomial subclass, we also perform a numerical experiment with

$$f(u) = |u|^{\zeta}, \quad g(v) = |v|^{\eta}. \qquad (4.1)$$

It turns out that, for either $\eta \le 2$ (irrespective of $\zeta$) or $\zeta \le 2$ (irrespective of $\eta$), the potential exhibits flat directions. This seems to allow for non-catastrophic runaway behaviour, the
<figure><img src="image_15.png" /><figcaption>**Figure 12.** We plot the ratio of the average stable evolution time $t_{\text{mean}}$ and the maximal chosen evolution time of $t_{\text{max}} = 10^2$ as a function of the polynomial degree $N$ in the polynomial subclass of integrable Liouville models, cf. eq. (3.68). For the coefficients, we choose $C_N = 1$ and $C_{n<N} = 0$ such that all even cases are stable according to section 3.3. Each point in the plot represents an average of the evolution of $N_{\text{evol}} = 100$ random initial conditions drawn from a normal distribution with standard deviation SD = 1.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** We plot the ratio of the average stable evolution time $t_{\text{mean}}$ and the maximal chosen evolution time of $t_{\text{max}} = 10^3$ as a function of (i) the exponent $\eta$ (at fixed $\zeta = 2$) in the left-hand panel, and for (ii) the exponent $\zeta$ (at fixed $\eta = 2$) in the right-hand panel, cf. eq. (4.1). Each point in the plot represents an average of the evolution of $N_{\text{evol}} = 100$ random initial conditions drawn from a normal distribution with standard deviation SD = 1. The numerical evolution is terminated if stiffness is detected or if $x(t) > C$ or $y(t) > C$, with $C = 10^3$.</figcaption></figure>

detection of which requires criterion (B). Without criterion (B), we do not detect any runaway behaviour for $\eta \le 2$ or $\zeta \le 2$ which indicates that the runaways occur with polynomial growth rate. We vary $\eta$ at fixed $\zeta = 2$ and vice versa: the results are shown in figure 13 and reproduce the conditions required in the proof of section 3.1. We note that close to the critical case, the potential becomes almost flat and the runaways are, therefore, increasingly hard to confidently detect numerically.
## 4.3 Numerical evidence in the absence of integrability

Going beyond integrability, we specify the following two classes of interaction potentials, i.e.,

$$ \tilde{V}_{I,\delta} = (1 + x^2 + y^2)^{1-\delta}, \qquad (4.2) $$

$$ \tilde{V}_{II,\delta} = (1 + x^2 y^2)^{\frac{1}{2}-\delta}, \qquad (4.3) $$

both supplemented with standard quadratic terms such that

$$ V_{I,\delta} = \tilde{V}_{I,\delta} + \frac{1}{2}x^2 - \frac{1}{2}y^2, \qquad (4.4) $$

$$ V_{II,\delta} = \tilde{V}_{II,\delta} + \frac{1}{2}x^2 - \frac{1}{2}y^2. \qquad (4.5) $$

(For simplicity, we only consider equal frequencies for $x$ and $y$.) In both cases, the exponent $\delta$ is chosen such that, as expected from the statement (N), $\delta \le 0$ determines the presence/absence of runaways. This is also what we find numerically (see below).

The above two example potentials involve non-polynomial interactions (except for special choices of $\delta$). Our third choice of example potential starts with a strictly polynomial interaction, i.e.,

$$ V_{\text{int}}(x, y) \equiv \tilde{V}_{III,\delta} = x^2 y^2. \qquad (4.6) $$

Supplemented only with standard kinetic terms, this interaction is unstable. However, we can generalize the decoupled potentials to

$$ V_P = (x^2)^{2+\delta}, \qquad V_G = -(y^2)^{2+\delta}, \qquad (4.7) $$

in which case, we find that for the full potential, i.e.,

$$ V_{III,\delta} = \tilde{V}_{III,\delta} + (x^2)^{2+\delta} - (y^2)^{2+\delta}, \qquad (4.8) $$

the exponent $\delta \le 0$, once more, seems to determine the presence/absence of runaways (at least within our finite evolution time).

We show the result of the respective numerical experiments in figure 14 for $V = V_{I,\delta}$ (upper panel), for $V = V_{II,\delta}$ (middle panel), and for $V = V_{III,\delta}$ (lower panel). In all of these cases, the inclusion of criterion (B) does not affect the results. The latter suggests that, for $\delta < 0$, the runaways are catastrophic, i.e., occur with exponentially fast growth rate.

Model III demonstrates in what sense added self-interactions can apparently stabilize the motion of nonintegrable systems or, at least, dramatically increase the longevity of stable motion. In the following, we give further detail on the apparent absence of runaways and the respective characterization of the motion as longlived.

We start by noting that the numerical experiment summarized in the lower panel of figure 14 shows that for $\delta < 0$, a significant fraction of initial conditions leads to a catastrophic runaway before $t_{\text{max}} = 100$.

This is in stark contrast to all cases with $\delta > 0$, for which we observe no sign of instabilities at all. The specific case of $\delta = 1$ is further exemplified in figure 15 where we evolve two
<figure><img src="image_17.png" /><figcaption>**Figure 14.** We plot the ratio of the average stable evolution time $t_{\text{mean}}$ and the maximal chosen evolution time of $t_{\text{max}} = 100$ as a function of the power-law exponent $\delta$ for the three classes of potentials $V_{I,\delta}$ (upper panel, cf. eq. (4.4)), $V_{II,\delta}$ (middle panel, cf. eq. (4.5)), and $V_{III,\delta}$ (lower panel, cf. eq. (4.8)). Each point in the plot represents an average of the evolution of $N_{\text{evol}} = 100$ random initial conditions drawn from a normal distribution with standard deviation SD = 1.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 15.** We show the evolution of two exemplary cases of random initial conditions (drawn from a Gaussian distribution with SD = 1 but selected with regards to the discussion in the main text) for model III with $\delta = 1$. In the upper panels, the motion remains in the local region in which non-integrable interactions dominate: hence the motion is apparently chaotic on all timescales. In the lower panels, the motion transitions to larger phase-space variables for which the decoupled — and thus approximately integrable — potential dominates: hence the motion exhibits seemingly non-chaotic phases (cf. lower left-hand panel). At very long timescales, both initial conditions lead to chaotic motion.</figcaption></figure>

distinct and random initial conditions up to $t_{max} = 4 \times 10^4$. We generically observe two different types of evolution, depending on the specified initial conditions.

On the one hand, initial conditions which are sufficiently small, lead to motion which remains within the central region, throughout the entire evolution. In these cases, the non-integrable character is always apparent. An example is shown in the upper panels of figure 15.

On the other hand, initial conditions which are sufficiently large, lead to motion which is dominated by the decoupled — and thus approximately integrable — self-interactions. In these cases, the motion is close to integrable. In general we find that, the larger the initial conditions, the harder it is to reveal the nonintegrable character of the motion. We refrain from explicitly showing one of the cases for which we cannot resolve the nonintegrable character up to $t_{max} = 4 \times 10^4$.

We also observe that, for some intermediate cases of initial conditions, the motion can transition between both types of behaviour. An example of this is shown in the lower panels of figure 15. At apparently random evolution times, the motion transitions back and forth between the nonintegrable character in the central region and the approximately integrable character at larger phase-space values.
We emphasize once more that, in all of these cases, the motion remains bounded in a
finite phase-space region, at least up to the investigated evolution time.

5 Discussion

Two point particles with opposite-sign kinetic energies can interact and yet exhibit stable motion. We prove this statement for a large class of integrable models, for which we show boundedness of motion, irrespective of the initial conditions (Lagrange stability). We gather substantial numerical evidence that nonintegrable models can be stable as well or at least longlived: within the necessarily finite numerical evolution time, we detect bounded motion and no indication of growth modes. We observe that our integrable and nonintegrable results seem to be connected by physical conditions on the interaction potential, which can be summarized as follows (as expressed in (N) of section 4): suppose a coupled Hamiltonian point-particle system with opposite-sign kinetic terms can be decoupled, in such a way that uncoupled ghosts and positive energy degrees of freedom are both separately stable due to their own self-interactions. Then, the coupled system avoids runaways, or at least exhibits very long-lived motion, if the full potential at sufficiently large coordinates is sufficiently dominated everywhere by its decoupled counterpart.

Our analytical results are contained in sections 2 and 3 and concern integrable two-particle Hamiltonians with an additional first integral which is quadratic in the momenta. (All the corresponding known cases with positive energies are collected in appendix A.) The analytical results are organized as follows.

First, we show how one can relate ghostly and non-ghostly theories by way of complex (univalent) canonical transformations and use this to build an integrable ghostly theory from a one with no ghost, i.e. to “ghostify” a positive kinetic energy theory while maintaining its integrability. We note that such ghostifications are *not* a mere recasting of variables as we do *not* complexify the initial conditions accordingly. Rather, real initial conditions on one side of the pair of theories correspond to purely complex initial conditions on the other side, and vice versa, while we work only with real variables when considering a given theory in isolation. Nevertheless, the ghostification procedure preserves integrability and thus serves as a generic construction principle for integrable models with ghostly interactions. Depending on the specific case, the Hamiltonian on either side of the pair of ghostly and non-ghostly theories can be complex-valued (and thus of little physical interest) while the other side may or may not be real-valued. With regards to the systematic construction of integrable ghostly Hamiltonians, two cases are of particular interest: first, the case in which both sides of the pair are real-valued. Second, the case in which a real-valued ghostly Hamiltonian can be obtained from a complex-valued non-ghostly Hamiltonian. We discuss examples of both these cases.

Second, we specify to a class of models with two free functions of (specific combinations) of the position variables, first obtained by Liouville. We delineate conditions on the free functions of the Liouville class and — assuming that said conditions hold — provide a rigorous proof that the resulting motion is bounded in a finite region in phase space, irrespective of the initial conditions. As in [16], the key to the proof is the additional integral of motion which, in combination with the Hamiltonian, can be used to bound the absolute value of all phase-space variables. This integral of motion is constructed from a positive definite quadratic
form of the canonical momenta, similar to the usual kinetic energy, and a “potential” part which is bounded from below and is unbounded from above. The conditions enunciated on the free functions ensure a sufficiently fast growth rate of the “potential” at large phase-space variables. These properties are in direct analogy with the usual positive definite conserved energy usually employed to establish stability. Note that for this question it is irrelevant which of the two conserved quantities drives the evolution. We thus expect that other integrable models can be analysed in a similar manner.

Third, we systematically identify an exhaustive polynomial subclass of the Liouville model. This polynomial subclass is real-valued both for the ghostly and the non-ghostly case. To the best of our knowledge, it has not been previously identified in the literature. As for stable motion with ghostly interactions, this polynomial subclass contains, in particular, a tower of polynomial potentials for which the proof conditions hold and for which we have thus established boundedness of motion in the presence of a ghost. The lowest non-trivial order features a sextic interaction term and quadratic terms which can be chosen at will. This provides a strikingly simple model of ghostly interactions which is polynomial, integrable, and globally stable. Moreover, this system exhibits spontaneous symmetry breaking entirely due to the interactions with the ghost. Remarkably, the very same interactions which are often thought to catastrophically destabilize the system, instead introduce novel (Lyapunov) stable vacua.

We highlight that all the above three parts of our analytical work are systematically presented and can, presumably, be extended to other integrable systems. Our analysis is therefore far from exhaustive: even among the known set of all integrable two-particle systems with a constant of motion, which is quadratic in the momenta (see appendix A), it is thus possible that further stable ghostly systems can be obtained in the same way.

Our numerical results are contained in section 4 and concern the extension to nonintegrable models. The analytical proof for boundedness of motion in the integrable Liouville class is quite striking. It is, however, equally striking that we find no indication that integrability is necessary. All of the investigated nonintegrable models seem to behave like the integrable ones. If the interactions at large phase-space variables do not obey the observed conditions (¤), we quickly detect runaway solutions which either lead to a breakdown of the evolution in finite time or to a less severe runaway. In contrast, whenever the interactions at large phase-space variables obey the observed conditions (¤), we find no indication for the onset of runaways.

The present work calls for various extensions in obvious directions. First, it is clear that our results should allow for an extension to the multiparticle case, i.e., to one or several ghost DoFs interacting with several positive-kinetic-energy DoFs. Second, an obvious question is that of quantization. A similar issue has been studied for some specific analogous PG system where stability could only be studied numerically and/or where the ghost still caused an instability, albeit an arguably “benign”, i.e., a slowly-growing one (see e.g. [29]). In this context the simple sextic polynomial potential which we have established to be integrable and stable presents a starting point to investigate whether stability persists, as we actually expect, at the quantum level. Last, we are planning to extend our work to field theories. In particular, one can construct a theory of two interacting scalar fields $\phi$ and $\psi$, one with positive and one with negative kinetic energy, for which the interaction potential is equivalent
to the one investigated here. Once more, the simple sextic polynomial potential found here, appears to be especially amenable for such a study. Finally, an extension of our analysis to gravitational interactions remains an important question for future work.

## Acknowledgments

AH would like to thank Claudia de Rham for discussion. CD thanks Thibault Damour and Robert Conte for discussions. We would like to thank Antonios Mitsopoulos for his helpful correspondence regarding appendix A. AH acknowledges: the work leading to this publication was supported by the PRIME programme of the German Academic Exchange Service (DAAD) with funds from the German Federal Ministry of Education and Research (BMBF). AH also acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) under Grant No 406116891 within the Research Training Group RTG 2522/1. AV is supported in part by the European Regional Development Fund (ESIF/ERDF) and the Czech Ministry of Education, Youth and Sports (MŠMT) through the Project CoGraDS — CZ.02.1.01/0.0/0.0/15 003/0000437. SM and AV collaboration is also supported by the Bilateral Czech-Japanese Mobility Plus Project JSPS-21-12 (JPJSBP120212502).

## A Integrable two-particle Hamiltonians with polynomial constants of motion up to quadratic order in the momenta

Here, we collect all the known pairs of a potential $V(x, y)$ and the respective constant of motion $I(x, y, p_x, p_y)$ which correspond to integrable Hamiltonian two-particle systems (cf. eq. (2.6)) and for which $I$ is either linear or quadratic in the momenta. A systematic search for these integrable models has been pioneered [40], formalized in [41], and subsequently completed in [42–49]. All results stated here can also be found throughout [49]. For the convenience of the reader (and correcting some typos), we collect them here in compact form. In order to match with the literature, we provide the integrable classes for the PP case. Using the “ghostification” procedure outlined in section 2, one can obtain the corresponding PG classes.

At linear order, only one class of integrable systems is found, i.e.,

$$ \begin{array}{l} \text{class 0: } V = f\left(\frac{a}{2}(x^2 + y^2) + cx - by\right), \qquad (A.1) \\ I = a(yp_x - xp_y) - bp_x - cp_y, \end{array} $$

with $a, b, c \in \mathbb{C}$ and $f$ an arbitrary function.

At quadratic order, we follow [49] and introduce shorthand variables (repeating some of the definitions given in the main text)

$$ r^2 = x^2 + y^2, \qquad (A.2) $$

$$ z_{\pm} = x \pm iy, \qquad (A.3) $$

$$ u^2 = \frac{1}{2} \left( r^2 + c + \sqrt{(r^2 + c)^2 - 4cx^2} \right), \qquad (A.4) $$
$$v^2 = \frac{1}{2} \left( r^2 + c - \sqrt{(r^2 + c)^2 - 4cx^2} \right), \qquad (A.5)$$

$$\hat{u}^2 = \frac{1}{2} \left( r^2 + \sqrt{r^4 - 2cz_{\pm}^2} \right), \qquad (A.6)$$

$$\hat{v}^2 = \frac{1}{2} \left( r^2 - \sqrt{r^4 - 2cz_{\pm}^2} \right). \qquad (A.7)$$

With these definitions the eight distinct classes of integrable systems can be written as

$$\text{class 1: } V = \frac{f(u) - g(v)}{u^2 - v^2}, \qquad (A.8)$$

$$I = -(xp_y - yp_x)^2 - c p_x^2 + 2 \frac{u^2 g(v) - v^2 f(u)}{u^2 - v^2},$$

$$\text{class 2: } V = g(r) + \frac{f(x/y)}{r^2}, \qquad (A.9)$$

$$I = (xp_y - yp_x)^2 + 2 f(x/y),$$

$$\text{class 3: } V = \frac{f(\hat{u}) - g(\hat{v})}{\hat{u}^2 - \hat{v}^2}, \qquad (A.10)$$

$$I = -(xp_y - yp_x)^2 - \frac{c}{2} (p_x \pm i p_y)^2 + 2 \frac{\hat{u}^2 g(\hat{v}) - \hat{v}^2 f(\hat{u})}{\hat{u}^2 - \hat{v}^2},$$

$$\text{class 4: } V = \frac{f(r + y) + g(r - y)}{r}, \qquad (A.11)$$

$$I = (yp_x - xp_y)p_x + \frac{(r+y)g(r-y) - (r-y)f(r+y)}{r},$$

$$\text{class 5: } V = \frac{f(z_+ + \sqrt{z_-}) + g(z_+ - \sqrt{z_-})}{\sqrt{z_-}}, \qquad (A.12)$$

$$I = (yp_x - xp_y)(p_x + i p_y) + i \frac{(p_x - i p_y)^2}{8} - i \left( \frac{z_+ - \sqrt{z_-}}{\sqrt{z_-}} \right) f(z_+ + \sqrt{z_-}) - i \left( \frac{z_+ + \sqrt{z_-}}{\sqrt{z_-}} \right) g(z_+ - \sqrt{z_-}),$$

$$\text{class 6: } V = \frac{f(z_{\pm})}{r} + g'(z_{\pm}), \qquad (A.13)$$

$$I = (yp_x - xp_y)(p_x \pm i p_y) \pm ig(z_{\pm}) \mp i z_{\pm} \left( \frac{f(z_{\pm})}{r} + g'(z_{\pm}) \right),$$

$$\text{class 7: } V = f(x) + g(y), \qquad (A.14)$$

$$I = p_x^2 + 2f(x),$$

$$\text{class 8: } V = r^2 f''(z_{\pm}) + g(z_{\pm}), \qquad (A.15)$$

$$I = p_x(p_x \pm i p_y) + r^2 f''(z_{\pm}) + g(z_{\pm}) + 2(z_{\pm})f'(z_{\pm}) - 2f(z_{\pm}),$$

where $c \in \mathbb{C}$ and $f$ and $g$ denote arbitrary functions. A prime denotes a first derivative and a double prime denotes a second derivative with respect to the argument. In comparison to [49], we have corrected a typo in the integral of motion $I$ of class 3 (given here in eq. (A.10)) which differs from [49, eq. (3.2.13)] but agrees with [55, eq. (41)] (given there in terms of $A = c/2$, $F_1 = 2f$, and $F_2 = 2g$). Invariants which include linear and quadratic terms in the momenta are prohibited by time reflection symmetry, cf. [49].
The above functional classes are known to contain several polynomial subclasses, cf. [49] for a collection of those. In appendices D and E, we identify the polynomial subclasses for class 1 (in the main text also referred to as Liouville class) and class 3. To the best of our knowledge these polynomial subclasses have not been previously identified.

## B Integrability and vacua

For self-completeness of the paper, as well as to obtain some intermediate results useful for the discussion of section 3.5, let us briefly review when a system with Hamiltonian

$$H = \frac{1}{2} p_x^2 + \frac{\sigma}{2} p_y^2 + V(x, y), \qquad (B.1)$$

possess a first integral that is quadratic in momenta

$$I = K^{ik}(x, y)p_i p_k + U(x, y), \qquad (B.2)$$

where $i, k$ run over $\{x, y\}$ and $\sigma = \pm 1$. The vanishing of the Poisson bracket $\{H, I\}_{\text{PB}}$ implies that

$$\begin{aligned} p_x[\partial_x U - 2(K^{xx} \partial_x V + K^{xy} \partial_y V)] + & \qquad (B.3) \\ & + p_y[\sigma \partial_y U - 2(K^{yx} \partial_x V + K^{yy} \partial_y V)] + \\ & + p_i p_k (p_x \partial_x K^{ik} + \sigma p_y \partial_y K^{ik}) = 0. \end{aligned}$$

This equation can only be satisfied for arbitrary $p_x$ and $p_y$ provided the following conditions hold

$$\begin{aligned} \partial_x U &= 2(K^{xx} \partial_x V + K^{xy} \partial_y V), \qquad (B.4) \\ \partial_y U &= 2\sigma(K^{yx} \partial_x V + K^{yy} \partial_y V), \\ p_i p_k (p_x \partial_x K^{ik} + \sigma p_y \partial_y K^{ik}) &= 0. \end{aligned}$$

As we assume that $K^{ik}(x, y)$ is a smooth function of $x$ and $y$, the first two equations of this system in eq. (B.4) imply that extrema of $V(x, y)$ are also extrema of $U(x, y)$. In particular, this means that equilibrium points of the system in eq. (B.1) are always located at the extrema of $U$. The last equality of the system in eq. (B.4) can be written as

$$p_x^3 \partial_x K^{xx} + p_x^2 p_y (\sigma \partial_y K^{xx} + 2 \partial_x K^{xy}) + p_y^2 p_x (\partial_x K^{yy} + 2 \sigma \partial_y K^{xy}) + \sigma p_y^3 \partial_y K^{yy} = 0. \qquad (B.5)$$

Again, this equation should hold for arbitrary $p_x$ and $p_y$, which results in the following conditions

$$\partial_x K^{xx} = \partial_y K^{yy} = 0, \qquad (B.6)$$

$$\partial_y K^{xx} = -2\sigma \partial_x K^{xy}, \quad \text{and} \quad \partial_x K^{yy} = -2\sigma \partial_y K^{xy}, \qquad (B.7)$$

where we used that $\sigma^2 = 1$ to go from eq. (B.5) to eq. (B.7). Differentiating the first of the conditions in eq. (B.7) with $x$ and second with $y$ and using eq. (B.6) we obtain

$$\partial_x^2 K^{xy} = \partial_y^2 K^{xy} = 0, \qquad (B.8)$$
which yields

$$K^{xy} = \kappa_0 + \kappa_x x + \kappa_y y + \kappa_{xy} xy, \qquad (B.9)$$

where $\kappa_0$, $\kappa_x$, $\kappa_y$ and $\kappa_{xy}$ are constants. Substituting this expression into the conditions in eq. (B.7) provides

$$\partial_y K^{xx} = -2\sigma(\kappa_x + \kappa_{xy}y),$$

and consequently (using also eq. (B.6)).

$$K^{xx} = \kappa_{xx} - 2\sigma \left( \kappa_x y + \frac{1}{2} \kappa_{xy} y^2 \right), \qquad (B.10)$$

and analogously

$$K^{yy} = \kappa_{yy} - 2\sigma \left( \kappa_y x + \frac{1}{2} \kappa_{xy} x^2 \right), \qquad (B.11)$$

where $\kappa_{xx}$ and $\kappa_{yy}$ are constants. Hence, the kinetic part can be also written as

$$K^{ik} p_k p_i = 2[\kappa_0 + \kappa_x x + \kappa_y y]p_x p_y + [\kappa_{xx} - 2\sigma\kappa_x y]p_x^2 + [\kappa_{yy} - 2\sigma\kappa_y x]p_y^2 - \sigma\kappa_{xy}[yp_x - \sigma xp_y]^2. \qquad (B.12)$$

The integral of motion $I$ is useful to establish stability when the kinetic part of $I$ is positive definite. This can be achieved when the terms linear in coordinates are absent i.e. when $\kappa_x = \kappa_y = 0$. Let us compare the expression in eq. (B.12) with the class 1 integral of motion $J_{LV}$ from eq. (3.31) with $\sigma = -1$ where

$$K_{LV}^{ik} p_k p_i = (xp_y + yp_x)^2 + \frac{\tilde{c}}{2} (p_x^2 + p_y^2),$$

we have $\kappa_0 = \kappa_x = \kappa_y = 0$ and the non-vanishing coefficients are

$$\kappa_{xx} = \kappa_{yy} = \frac{\tilde{c}}{2}, \quad \text{and} \quad \kappa_{xy} = 1.$$

After the form of $K^{ik}$ is fixed, one can solve the system of first-order PDEs in eq. (B.4) to find $V(x, y)$ and $U(x, y)$. This task is beyond the scope of the paper. The corresponding solutions of this system are reviewed in e.g [49] see also appendix A.

The main point of this appendix was to demonstrate that equilibrium configurations $(\partial_x V = \partial_y V = 0)$ of integrable systems as in eq. (B.1) are all located at extrema of $U(x, y)$. However, the opposite is not true in general. Indeed, as the structure of $K^{ik}$ in eq. (B.12) clearly shows, the matrix of the system of the first two equations in eq. (B.4) may not be invertible. In particular, one can choose $\kappa_{xy} = \kappa_y = \kappa_{yy}$ so that $K^{yy} = 0$ and one cannot invert eq. (B.4) to uniquely express $\partial_i V$ through $\partial_k U$. This happens, e.g., for the models in eqs. (A.11) and (A.15).

Finally, let us discuss a useful relation between the Hessian determinant of the potential $U$,

$$\mathcal{H}^{\text{ess}} U = \partial_x^2 U \partial_y^2 U - \partial_x \partial_y U \partial_x U, \qquad (B.13)$$
and that of the potential $V$ at an equilibrium point. The vanishing of the derivatives of $V$ at such points implies, using eq. (B.4), that the derivatives of $K^{ik}$ will make no contribution to second derivatives of $U$ there. Using then the first equation in (B.4) to compute $\partial_y \partial_x U$ and the second equation in (B.4) to compute $\partial_x \partial_y U$, which facilitates the computation at equilibrium points, one finds after some algebra, that one has there

$$ \mathcal{H}_{\text{ess}} U = 4\sigma (\det K) \mathcal{H}_{\text{ess}} V. \qquad (B.14) $$

If the first integral is useful to establish stability, then $K^{ik}$ is positive definite and not degenerate, so that $\det K > 0$. In that case, for PP case $\sigma = 1$ so that the minimum of $U$ corresponds to the minimum of $V$. However, for PG case $\sigma = -1$ and $\mathcal{H}_{\text{ess}} U > 0$ corresponds to $\mathcal{H}_{\text{ess}} V < 0$. Hence the minimum of $U$ can only be achieved at saddle points of $V$. Thus, stable equilibrium points in systems where a ghost $y$ interacts with a usual DoF $x$ are located in saddle points of the potential $V$.

### C Stable motion with a ghostly $x$

The proof of section 3.1 can easily be used to show that stable motions also exist in another class of PG Liouville models: the one obtained from the PP where one now ghostifies $x$ by a complex canonical transformation $\mathcal{E}_x^{\pm}$. I.e. starting from the PP Liouville model eq. (3.1) with $\sigma = +1$, we change $x$ to $ix$ and $p_x$ to $-ip_x$. This yields a model, dubbed here a GP model, with a Hamiltonian and an extra constant of motion respectively given by

$$ H_{\text{LV,GP}} = -\frac{p_x^2}{2} + \frac{p_y^2}{2} + V_{\text{LV}}(x, y) \qquad (C.1) $$

$$ I_{\text{LV,GP}} = (p_y x + p_x y)^2 + c p_x^2 + \mathcal{V} \qquad (C.2) $$

where $V_{\text{LV}}$, $\mathcal{V}$, $u$ and $v$ are defined respectively as in eqs. (3.2) to (3.4) and (3.7), and where $r^2$ now reads

$$ r^2 = -x^2 + y^2. \qquad (C.3) $$

Note in particular, that the above $I_{\text{LV}}$ is just identical to eq. (C.2) upon the replacement of the positive $\tilde{c}$ by a positive $c$. The analog of the constant of motion eq. (3.31) is now given by the linear combination $I_{\text{LV,GP}} + cH_{\text{LV,GP}}$. In fact, one can show the mere equivalence between the PG models (i.e. the one obtained from the Liouville PP model by applying $\mathcal{E}_y^{\pm}$) and the GP models (i.e. the one obtained from the Liouville PP model by applying $\mathcal{E}_x^{\pm}$). Indeed, considering the expressions for $u$ and $v$ with $\sigma = +1$ given in eqs. (3.3) and (3.4) and reminded below

$$ u^2 = \frac{1}{2} \left( r^2 + c + \sqrt{(r^2 + c)^2 - 4cx^2} \right), \qquad (C.4) $$

$$ v^2 = \frac{1}{2} \left( r^2 + c - \sqrt{(r^2 + c)^2 - 4cx^2} \right), \qquad (C.5) $$

$$ r^2 = x^2 + y^2, \qquad (C.6) $$
we see that the application of $\mathcal{C}_y^{\pm}$ to these expressions just changes $r^2 = x^2+y^2$ to $r^2 = x^2-y^2$ and hence transforms $u^2$ and $v^2$ to

$$ u^2 = \frac{1}{2} (x^2 - y^2 + c + \sqrt{(x^2 - y^2 + c)^2 - 4cx^2}), \qquad (C.7) \\ \equiv u_{\text{PG}}^2(c), $$

$$ v^2 = \frac{1}{2} (x^2 - y^2 + c - \sqrt{(x^2 - y^2 + c)^2 - 4cx^2}), \qquad (C.8) \\ \equiv v_{\text{PG}}^2(c). $$

Conversely, the application of $\mathcal{C}_x^{\pm}$ to eqs. (C.4) and (C.5) leads to

$$ u^2 = \frac{1}{2}(-x^2 + y^2 + c + \sqrt{(-x^2 + y^2 + c)^2 + 4cx^2}), \qquad (C.9) \\ \equiv u_{\text{GP}}^2(c) = -v_{\text{GP}}^2(-c), $$

$$ v^2 = \frac{1}{2}(-x^2 + y^2 + c - \sqrt{(-x^2 + y^2 + c)^2 + 4cx^2}), \qquad (C.10) \\ \equiv v_{\text{GP}}^2(c) = -u_{\text{GP}}^2(-c). $$

This implies that, starting from a given PP Hamiltonian of the form

$$ H_{\text{LV,PP}} = \frac{p_x^2}{2} + \frac{p_y^2}{2} + \frac{\Phi(u^2) - \Psi(v^2)}{u^2 - v^2}, \qquad (C.11) $$

where for simplicity, we have assumed the functions $f$ and $g$ of eq. (3.2) to be respectively function $\Phi$ and $\Psi$ of $u^2$ and $v^2$, i.e. $f(u) = \Phi(u^2)$ and $g(v) = \Psi(v^2)$, we get by applying $\mathcal{C}_y^{\pm}$ the PG theory of Hamiltonian

$$ H_{\text{LV,PG}} = \frac{p_x^2}{2} - \frac{p_y^2}{2} + \frac{\Phi(u_{\text{PG}}^2(c)) - \Psi(v_{\text{PG}}^2(c))}{u_{\text{PG}}^2(c) - v_{\text{PG}}^2(c)}. \qquad (C.12) \\ \equiv H_{\text{LV,PG}}(c, \Phi, \Psi) $$

where $H_{\text{LV,PG}}$ is a function of $c$ but a functional of $\Phi : u^2 \mapsto \Phi(u^2)$ and $\Psi : u^2 \mapsto \Psi(u^2)$. On the other hand, the application of $\mathcal{C}_x^{\pm}$ on the same PP theory yields an Hamiltonian

$$ H_{\text{LV,GP}} = -\frac{p_x^2}{2} + \frac{p_y^2}{2} + \frac{\Phi(u_{\text{GP}}^2(c)) - \Psi(v_{\text{GP}}^2(c))}{u_{\text{GP}}^2(c) - v_{\text{GP}}^2(c)}, \qquad (C.13) \\ \equiv H_{\text{LV,GP}}(c, \Phi, \Psi), \\ = -\frac{p_x^2}{2} + \frac{p_y^2}{2} + \frac{-\Psi(-u_{\text{PG}}^2(-c)) + \Phi(-v_{\text{PG}}^2(-c))}{u_{\text{PG}}^2(-c) - v_{\text{PG}}^2(-c)}, $$

where, as above, $H_{\text{LV,GP}}$ is a function of $c$ but a functional of $\Phi : u^2 \mapsto \Phi(u^2)$ and $\Psi : u^2 \mapsto \Psi(u^2)$. So we have that

$$ H_{\text{LV,GP}}(c, \Phi, \Psi) = -H_{\text{LV,PG}}(-c, u^2 \mapsto \Psi(u^2), v^2 \mapsto \Phi(v^2)) \qquad (C.14) $$

showing that the PG and GP theories just differ by a global sign in the Hamiltonian, a proper redefinition of the function $\Phi$ and $\Psi$ (or $f$ and $g$), and a change of sign of $c$; hence these theories are equivalent. In particular, the conditions enunciated in (i) and (ii) on the functions $f$ and $g$ of the proof of section 3.1 yield a stable theory PG eq. (C.12) for $c < 0$ and a stable theory GP eq. (C.12) for $c > 0$.
## D Polynomial Liouville potentials

We have systematically constructed all potentials which are polynomials in $x$ and $y$ of degree $2N$ (with $N \in \mathbb{N}$) and belong to the class of integrable Liouville potentials. (In the following, we present them for the PG case, which is of relevance for this paper. With reference to section 2, these can straightforwardly be transformed to the PP case.)

At fixed degree $2N$, the polynomial potential can be found by (i) making an ansatz with arbitrary coefficients, (ii) transforming from $(x, y)$ to $(u, v)$ coordinates, (iii) multiplying by $u^2 - v^2$, and (iv) demanding that all mixed coefficients vanish, i.e., that the potential is of the integrable form in eq. (3.2). We solve this procedure up to order $N = 8$ and infer that the general form of these potentials corresponds to choosing

$$f(u) = \sum_{n=1}^{N} C_n (u^2)^n, \qquad (D.1)$$

$$g(v) = \sum_{n=1}^{N} C_n (v^2)^n. \qquad (D.2)$$

Proof that this leads to polynomial $V(x, y)$ is given as follows. We first note that $u^2 = \frac{1}{2}(r^2 + c + \tilde{W})$ while $v^2 = \frac{1}{2}(r^2 + c - \tilde{W})$ where $\tilde{W}$ is defined in eq. (3.16) With this definition, we may rewrite the potential in eq. (3.68) as follows

$$ \begin{aligned} V_{\text{LV}}^{(N)} &= \sum_{n=1}^{N} \frac{C_n}{\tilde{W}} \left[ (u^2)^n - (v^2)^n \right] && (\text{D.3}) \\ &= \sum_{n=1}^{N} \frac{C_n}{2^n \tilde{W}} \left[ \left( (r^2 + c) + \tilde{W} \right)^n - \left( (r^2 + c) - \tilde{W} \right)^n \right] \\ &= \sum_{n=1}^{N} \sum_{k=1}^{n} \frac{C_n}{2^n} \left[ \binom{n}{k} (r^2 + c)^{n-k} \tilde{W}^{k-1} (1 - (-1)^k) \right], \end{aligned} $$

where we have expanded and collected the binomial coefficients in the last step. Each summand with even $k$ vanishes due to the last factor. Each summand with odd $k$ gives a contribution in which $\tilde{W}$ is raised to an even power. Hence, all of the summands are indeed polynomials in $x$ and $y$.

Given that all summands with even $k$ vanish, we may also shift the sum: in terms of the variables $(x, y)$ the polynomial potential thus reads

$$V_{\text{LV}}^{(N)}(x, y) = \sum_{n=1}^{N} \sum_{k=1}^{\lfloor \frac{n}{2} \rfloor - 1} \left[ \frac{C_n}{2^{n-1}} \binom{n}{2k+1} (r(x,y)^2 + c)^{n-(2k+1)} (\tilde{W}(x,y)^2)^k \right], \quad (\text{D.4})$$

where the $\lceil \rceil$ designate the ceiling function which returns the lowest integer larger or equal to its argument. The explicit first few terms of the potential read

$$ \begin{aligned} V_{\text{LV}}^{(N)}(x, y) = &+ C_1 + C_2(x^2 - y^2 + c) + C_3 \left[ (x^2 - y^2 + c)^2 - c x^2 \right] \\ &+ C_4 \left[ (x^2 - y^2 + c)^3 - 2 c x^2 (x^2 - y^2 + c) \right] \\ &+ \dots, \end{aligned} \qquad (\text{D.5}) $$
from which we can see that the $n = 1$ term corresponds to a constant shift of the potential, the $n = 2$ term reproduces equal (but opposite-sign) quadratic terms for the $x$ and the $y$ mode, and higher-order $n \ge 3$ terms correspond to specific polynomials.

Generally, the $n^{th}$-order polynomial also includes (some of) the terms appearing at lower order. In particular, all $C_n$ contribute to the quadratic term for the $x$ mode. At fixed order, the quadratic terms can thus be adjusted freely. For instance, at order $n = 4$, we may trade $C_2$ and $C_3$ for the frequencies $\omega_x^2$ and $\omega_y^2$, as well as $C_1$ to remove the constant shift. More specifically, we define

$$ 
\begin{aligned}
C_1 &= -c \frac{\omega_x^2}{2}, \\
C_2 &= \omega_x^2 - \frac{\omega_y^2}{2} + c^2 C_4, \\
C_3 &= -\frac{\omega_x^2 - \omega_y^2}{2c} - 2c C_4.
\end{aligned} 
\qquad (D.6) 
$$

The $4^{th}$-order potential $V_4(x, y)$ can then be rewritten as

$$ V_{LV}^{(4)}(x, y) = \frac{\omega_x^2}{2} \left[ x^2 - \frac{(x^2 - y^2)^2}{c} \right] - \frac{\omega_y^2}{2} \left[ y^2 - \frac{(x^2 - y^2)^2}{c} \right] + C_4 \left[ (x^2 - y^2)^3 - c(x^4 - y^4) \right], \qquad (D.7) $$

which can be rewritten as in eq. (3.69).

For $N = 6$ we find (reabsorbing an additional constant shift into a redefinition of $C_1$)

$$ 
\begin{aligned}
V_{LV}^{(6)}(x, y) = V_4(x, y) &+ C_5(x^2 - y^2)^4 + c C_5(x^2 - y^2)^2(x^2 - 4y^2) + c^2 C_5(3y^4 - 2x^4) \\
&+ C_6(x^2 - y^2)^5 + c C_6(x^2 - y^2)^3(x^2 - 5y^2) + c^2 C_6(x^2 - y^2)(x^4 - 8x^2y^2 + 10y^4) \\
&+ c^3 C_6(6y^4 - 3x^4)
\end{aligned} 
\qquad (D.8) 
$$

We note that the finite-order construction from which we started out suggests that eq. (D.4) is exhaustive and captures all possible polynomial potentials in the integrable Liouville model.

Finally, it is also instructive to note that, with the choice in eqs. (D.1) and (D.2) (leading to the polynomial potentials eq. (D.4) with some fixed $N$), the quantity $\mathcal{V}$ defined in eq. (3.7) and entering into the constant of motion $I_{LV}$ is also polynomial and reads

$$ \mathcal{V}^{(N)} = -2 u^2 v^2 \sum_{n=1}^{N-1} \frac{C_{n+1}}{\tilde{W}} \left[ (u^2)^n - (v^2)^n \right] \qquad (D.9) $$

where, considering e.g. eq. (D.3), the second term on the right-hand side above is just the same as $-2u^2v^2 \times V_{LV}^{(N-1)}$ where the $C_n$ are replaced by $C_{n+1}$. For example, for $N = 4$, we get, using eq. (3.17) as well as eq. (D.5),

$$ \mathcal{V}^{(4)} = -2 c x^2 \left[ C_2 + C_3(x^2 - y^2 + c) + C_4 \left( (x^2 - y^2 + c)^2 - c x^2 \right) \right], \qquad (D.10) $$
where the coefficients $C_1, C_2, C_3$ can be replaced as in eq. (D.6) to yield the expression corresponding to the form in eq. (D.7), i.e.,

$$ \mathcal{V}^{(4)} = -2c x^2 \left[ \frac{\omega_x^2}{2} - \frac{\omega_x^2 - \omega_y^2}{2c} (x^2 - y^2) + C_4 \left( (x^2 - y^2)^2 - c x^2 \right) \right]. \quad (D.11) $$

Clearly, this also implies that the alternative constant of motion $J_{LV}$ used in the proof and defined in eq. (3.31) is polynomial since $J_{LV} = I_{LV} + c H_{LV} = (xp_y + yp_x)^2 - \frac{c}{2}(p_x^2 + p_y^2) + \mathcal{V} + c V_{LV}$, or equivalently $\mathcal{U} = \mathcal{V} + c V_{LV}$. For completeness, we state the resulting explicit fourth-order form, i.e.,

$$ \mathcal{U}^{(4)} = \frac{\omega_x^2 - \omega_y^2}{2} (x^4 - y^4) - c \left[ \frac{\omega_x^2}{2} x^2 + \frac{\omega_y^2}{2} y^2 \right] - c C_4 \left[ (x^2 - y^2)^2 (x^2 + y^2) - c(x^4 + y^4) \right], \quad (D.12) $$

written again in terms of the frequencies and $C_4$.

## E Polynomial class-3 potentials

Class 3 in the previous appendix A provides an example for an integrable two-particle Hamiltonian which is complex-valued in the PP case but when the y-mode is ghostified (via the complex transformation $\mathcal{E}_y^{\pm}$, cf. section 2, i.e. taking $y \to \pm iy$ and $p_y \to \mp ip_y$) the resulting PG theory is real-valued. In fact, the complexification is real-valued for both signs in $z_{\pm}$ (as part of the model definition in eq. (A.3)) and thus, in fact, generates two complex classes. These are distinguished with the occurrence of “$\pm$” below.

The same choice for the functions $f(\hat{u})$ and $g(\hat{v})$ as in appendix D, eqs. (D.1) and (D.2), also leads to a polynomial subclass. (Note, however, that $\hat{u}$ and $\hat{v}$ for class 3 differ from $u$ and $v$ for the Liouville class / class 1.)

The polynomial form in terms of $(x, y)$ can be derived in complete analogy to eq. (D.4) and reads

$$ V_{(\text{class 3})}^{(N)}(x, y) = \sum_{n=1}^{N} \sum_{k=1}^{\lfloor \frac{n}{2} \rfloor - 1} \left[ \frac{C_n}{2^{n-1}} \binom{n}{2k+1} r^{2(n-(2k+1))} \hat{W}^{2k} \right], \quad (E.1) $$

where $\lceil \rceil$ denotes the ceiling function, $r = \sqrt{x^2 - y^2}$, and $\hat{W} = \sqrt{r^4 - 2c(x \pm y)^2}$ as in the definition of class 3 in appendix A (after complexification).

The explicit first few terms of this potential read

$$ \begin{aligned} V_{(\text{class 3})}^{(N)}(x, y) = &+ C_1 + C_2(x^2 - y^2) + C_3 \left[ (x^2 - y^2)^2 - \frac{1}{2}c(x \pm y)^2 \right] \\ &+ C_4 \left[ (x^2 - y^2)^3 - c(x^2 - y^2)(x \pm y)^2 \right] \\ &+ \dots \end{aligned} \quad (E.2) $$
Introducing the frequencies and removing a constant shift, i.e., fixing $C_1 = 0$, $C_2 = \frac{\omega_x^2+\omega_y^2}{4}$, and $C_3 = -\frac{\omega_x^2-\omega_y^2}{2c}$, the $N = 4$ case can be written as

$$V_{\text{(class 3)}}^{(4)}(x, y) = \frac{\omega_x^2}{2} x^2 - \frac{\omega_y^2}{2} y^2 - \frac{\omega_x^2 - \omega_y^2}{4c} \left[ 2(x^2 - y^2)^2 + c(x^2 + y^2) - c(x \pm y)^2 \right] \\ - C_4 c(x^2 - y^2)(x \pm y)^2 + C_4 (x^2 - y^2)^3. \qquad (\text{E.3})$$

In the main text, we specify to $N = 4$ and to equal quadratic terms, i.e., to $\omega_x = \omega_y \equiv \omega$.

## References

[1] M. Ostrogradsky, *Mémoires sur les équations différentielles, relatives au problème des isopérimètres*, *Mem. Acad. St. Petersburg* **6** (1850) 385 [INSPIRE].

[2] J.-P. Bruneton and G. Esposito-Farese, *Field-theoretical formulations of MOND-like gravity*, *Phys. Rev. D* **76** (2007) 124012 [Erratum ibid. **76** (2007) 129902] [arXiv:0705.4043] [INSPIRE].

[3] S.V. Ketov, G. Michiaki and T. Yumibayashi, *Quantizing with a higher time derivative*, *arXiv:1110.1155* [INSPIRE].

[4] R.P. Woodard, *Ostrogradsky's theorem on Hamiltonian instability*, *Scholarpedia* **10** (2015) 32243 [arXiv:1506.02210] [INSPIRE].

[5] A. Pais and G.E. Uhlenbeck, *On Field theories with nonlocalized action*, *Phys. Rev.* **79** (1950) 145 [INSPIRE].

[6] T.D. Lee and G.C. Wick, *Finite Theory of Quantum Electrodynamics*, *Phys. Rev. D* **2** (1970) 1033 [INSPIRE].

[7] K.S. Stelle, *Renormalization of Higher Derivative Quantum Gravity*, *Phys. Rev. D* **16** (1977) 953 [INSPIRE].

[8] B. Grinstein, D. O'Connell and M.B. Wise, *The Lee-Wick standard model*, *Phys. Rev. D* **77** (2008) 025012 [arXiv:0704.1845] [INSPIRE].

[9] A.D. Linde, *The Inflationary Universe*, *Rept. Prog. Phys.* **47** (1984) 925 [INSPIRE].

[10] A.D. Linde, *The Universe Multiplication and the Cosmological Constant Problem*, *Phys. Lett. B* **200** (1988) 272 [INSPIRE].

[11] D.E. Kaplan and R. Sundrum, *A Symmetry for the cosmological constant*, *JHEP* **07** (2006) 042 [hep-th/0505265] [INSPIRE].

[12] R. Brandenberger and P. Peter, *Bouncing Cosmologies: Progress and Problems*, *Found. Phys.* **47** (2017) 797 [arXiv:1603.05834] [INSPIRE].

[13] R.R. Caldwell, *A Phantom menace?*, *Phys. Lett. B* **545** (2002) 23 [astro-ph/9908168] [INSPIRE].

[14] E. Di Valentino et al., *In the realm of the Hubble tension — a review of solutions*, *Class. Quant. Grav.* **38** (2021) 153001 [arXiv:2103.01183] [INSPIRE].

[15] J.M. Cline, S. Jeon and G.D. Moore, *The Phantom menaced: Constraints on low-energy effective ghosts*, *Phys. Rev. D* **70** (2004) 043543 [hep-ph/0311312] [INSPIRE].

[16] C. Deffayet, S. Mukohyama and A. Vikman, *Ghosts without Runaway Instabilities*, *Phys. Rev. Lett.* **128** (2022) 041301 [arXiv:2108.06294] [INSPIRE].
[17] S.W. Hawking and T. Hertog, *Living with ghosts*, Phys. Rev. D 65 (2002) 103515 [hep-th/0107088] [INSPIRE].

[18] 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 [arXiv:0706.0207] [INSPIRE].

[19] J. Garriga and A. Vilenkin, *Living with ghosts in Lorentz invariant theories*, JCAP 01 (2013) 036 [arXiv:1202.1239] [INSPIRE].

[20] A. Salvio and A. Strumia, *Agravity*, JHEP 06 (2014) 080 [arXiv:1403.4226] [INSPIRE].

[21] A. Smilga, *Classical and quantum dynamics of higher-derivative systems*, Int. J. Mod. Phys. A 32 (2017) 1730025 [arXiv:1710.11538] [INSPIRE].

[22] D. Becker, C. Ripken and F. Saueressig, *On avoiding Ostrogradski instabilities within Asymptotic Safety*, JHEP 12 (2017) 121 [arXiv:1709.09098] [INSPIRE].

[23] D. Anselmi, *Fakeons And Lee-Wick Models*, JHEP 02 (2018) 141 [arXiv:1801.00915] [INSPIRE].

[24] J.F. Donoghue and G. Menezes, *Unitarity, stability and loops of unstable ghosts*, Phys. Rev. D 100 (2019) 105006 [arXiv:1908.02416] [INSPIRE].

[25] C. Gross, A. Strumia, D. Teresi and M. Zirilli, *Is negative kinetic energy metastable?*, Phys. Rev. D 103 (2021) 115025 [arXiv:2007.05541] [INSPIRE].

[26] J.F. Donoghue and G. Menezes, *Ostrogradsky instability can be overcome by quantum physics*, Phys. Rev. D 104 (2021) 045010 [arXiv:2105.00898] [INSPIRE].

[27] A. Platania, *Causality, unitarity and stability in quantum gravity: a non-perturbative perspective*, JHEP 09 (2022) 167 [arXiv:2206.04072] [INSPIRE].

[28] E. Pagani, G. Tecchiolli and S. Zerbini, *On the Problem of Stability for Higher Order Derivatives: Lagrangian Systems*, Lett. Math. Phys. 14 (1987) 311 [INSPIRE].

[29] A.V. Smilga, *Benign versus malicious ghosts in higher-derivative theories*, Nucl. Phys. B 706 (2005) 598 [hep-th/0407231] [INSPIRE].

[30] N. Boulanger, F. Buisseret, F. Dierick and O. White, *Higher-derivative harmonic oscillators: stability of classical dynamics and adiabatic invariants*, Eur. Phys. J. C 79 (2019) 60 [arXiv:1811.07733] [INSPIRE].

[31] T. Damour and A. Smilga, *Dynamical systems with benign ghosts*, Phys. Rev. D 105 (2022) 045018 [arXiv:2110.11175] [INSPIRE].

[32] S.M. Carroll, M. Hoffman and M. Trodden, *Can the dark energy equation-of-state parameter w be less than -1?*, Phys. Rev. D 68 (2003) 023509 [astro-ph/0301273] [INSPIRE].

[33] I.B. Ilhan and A. Kovner, *Some Comments on Ghosts and Unitarity: The Pais-Uhlenbeck Oscillator Revisited*, Phys. Rev. D 88 (2013) 044045 [arXiv:1301.4879] [INSPIRE].

[34] M. Pavšič, *Pais-Uhlenbeck oscillator and negative energies*, Int. J. Geom. Meth. Mod. Phys. 13 (2016) 1630015 [arXiv:1607.06589] [INSPIRE].

[35] M. Pavšič, *Stable Self-Interacting Pais-Uhlenbeck Oscillator*, Mod. Phys. Lett. A 28 (2013) 1350165 [arXiv:1302.5257] [INSPIRE].

[36] M. Pavšič, *On negative energies, strings, branes, and braneworlds: A review of novel approaches*, Int. J. Mod. Phys. A 35 (2020) 2030020 [arXiv:2012.04976] [INSPIRE].

[37] J.L. Salle and S. Lefschetz, *Stability By Liapunov's Direct Methods with Applications*, Academic Press (1961).
[38] J. Liouville, Note sur l'intégration des équations différentielles de la Dynamique, présentée au Bureau des Longitudes le 29 juin 1853, *J. Math. Pures Appl.* **20** (1855) 137.

[39] V.I. Arnold, *Mathematical methods of classical mechanics*, Vol. 60, Springer Science & Business Media (2013).

[40] G. Darboux, *Sur un problème de mécanique*, *Bull. Sci. Math. Astron.* **2** (1878) 433.

[41] E.T. Whittaker, *A treatise on the analytical dynamics of particles and rigid bodies*, CUPCambridge University Press Archive (1964).

[42] J. Friš, Y.A. Smorodinskii, M. Uhlír and P. Winternitz, *Symmetry groups in classical and quantum mechanics*, *Sov. J. Nucl. Phys.* **4** (1967) 444.

[43] C.R. Holt, *Construction of new integrable hamiltonians in two degrees of freedom*, *J. Math. Phys.* **23** (1982) 1037.

[44] A. Ankiewicz and C. Pask, *The complete whittaker theorem for two-dimensional integrable systems and its application*, *J. Phys. A* **16** (1983) 4203.

[45] B. Dorizzi, B. Grammaticos and A. Ramani, *A new class of integrable systems*, *J. Math. Phys.* **24** (1983) 2282.

[46] B. Grammaticos, B. Dorizzi and A. Ramani, *Integrability of Hamiltonians with third- and fourth-degree polynomial potentials*, *J. Math. Phys.* **24** (1983) 2289.

[47] G. Thompson, *Darboux's problem of quadratic integrals*, *J. Phys. A* **17** (1984) 985.

[48] T. Sen, *Integrable potentials with quadratic invariants*, *Phys. Lett. A* **111** (1985) 96.

[49] J. Hietarinta, *Direct Methods for the Search of the Second Invariant*, *Phys. Rept.* **147** (1987) 87 [inSPIRE].

[50] J. Liouville, *Sur quelques cas particuliers où les équations du mouvement d'un point matériel peuvent s'intégrer*, *J. Math. Pures Appl.* **11** (1846) 345.

[51] D. Robert and A.V. Smilga, *Supersymmetry vs ghosts*, *J. Math. Phys.* **49** (2008) 042104 [math-ph/0611023] [inSPIRE].

[52] W.R. Inc., *Mathematica, Version 13.1, ExplicitRungeKutta Method for NDSolve*, Champaign, IL, U.S.A. (2022) [https://reference.wolfram.com/language/tutorial/NDSolveExplicitRungeKutta.html].

[53] L.F. Shampine, *Numerical solution of ordinary differential equations*, Routledge (2018).

[54] W.R. Inc., *Mathematica, Version 13.1, Stiffness Detection*, Champaign, IL, U.S.A. (2022) [https://reference.wolfram.com/language/tutorial/NDSolveStiffnessTest.html].

[55] A. Mitsopoulos, M. Tsamparis and A. Paliathanasis, *Integrable and Superintegrable Potentials of 2d Autonomous Conservative Dynamical Systems*, *Symmetry* **12** (2020) 1655.