Post-publication activity

Curator: Richard P Woodard

Ostrogradsky's construction of a Hamiltonian formalism for nondegenerate higher derivative Lagrangians is reviewed. The resulting instability imposes by far the most powerful restriction on fundamental, interacting, continuum Lagrangian field theories. A discussion is given of the problems raised by attempts to evade this restriction.

Introduction

Albert Einstein famously commented, “What really interests me is whether God had any choice in the creation of the world” (Straus, 1956). Within the context of local Lagrangian field theory the answer seems to be that powerful restrictions exist but some freedom still remains regarding the choice of dynamical variables and symmetries. By far the greatest restriction is the obstacle to including higher time derivatives which is implied by Ostrogradsky's construction of a canonical formalism for nondegenerate higher derivative Lagrangians (Ostrogradsky, 1850).

Mikhail Vasilevich Ostrogradsky lived from 1801 to 1862. He was born to a poor family of Ukrainian ethnicity in Pashennaya, which is now in Ukraine but was at that time part of the vast Russian Empire. These were momentous years for Russia, bracketed by its rise to become the predominant military power during the Napoleonic Wars, and its humiliating collapse before Britain and France during the Crimean War. Russian society was riven by the struggle between the forces of reaction and reform. Indeed, Ostrogradsky was denied his doctorate at the University of Kharkov because the mathematics professor who had examined him was considered insufficiently religious. Later on, Ostrogradsky was placed under police surveillance at the start of his career in the Imperial Russian capital of St. Petersburg (O'Connor and Robertson).

Ostrogradsky studied and worked in Paris from 1822 through 1827. He knew the leading French mathematicians of the time, including Cauchy, who paid off his debts and secured him a teaching job. In 1826 Ostrogradsky stated and proved the divergence theorem, which was later re-discovered by Gauss in the 1830's. Ostrogradsky paid a much shorter visit to Paris in 1830. However, most of his professional life was spent in St. Petersburg where he was elected to the Imperial Academy of Sciences and played an important role in the teaching of mathematical sciences. Ostrogradsky wrote in French and Russian (O'Connor and Robertson).

Ostrogradsky's higher derivative generalization of Hamilton's construction was published in 1850 (Ostrogradsky, 1850). Ostrogradsky's construction implies that there is a linear instability in the Hamiltonians associated with Lagrangians which depend upon more than one time derivative in such a way that the higher derivatives cannot be eliminated by partial integration. This is probably why Newton was right to assume the laws of physics take the form of second differential equations when expressed in terms of fundamental dynamical variables.

It might seem curious that Ostrogradsky did not appreciate the importance of his construction to fundamental theory. However, one must recall that the researchers of his time were just beginning to make the connection between energy functionals and the concept of stability -- which in those days meant the absence of growing perturbations. The notion of quantum fluctuations exploring all perturbations was decades away, and the key insight that all dynamics is described by interacting continuum field theories was even further in the future.

The next section presents Ostrogradsky's construction in the context of point particle whose position is $x(t)$. The consequent sections discuss the consequences of this result for fundamental theory and deal with quantization and degeneracy. Lastly, we end with some concluding remarks.

This section presents Ostrogradsky's construction. First, the usual case of a first derivative Lagrangian is reviewed to fix concepts and notation. Then the case of second derivatives is presented. The section closes with a review of the general case of $N$ time derivatives.

Hamilton's Construction

In the usual case of $L = L(x,\dot{x})$, the Euler-Lagrange equation is, $\tag{1} \frac{\partial L}{\partial x} - \frac{d}{dt} \frac{\partial L}{\partial \dot{x}} = 0 \; .$ The assumption that $\frac{\partial^2 L}{\partial \dot{x}^2} \neq 0$ is known as nondegeneracy. If the Lagrangian is nondegenerate one can write (1) in the form Newton assumed so long ago for the laws of physics, $\ddot{x} = \mathcal{F}(x,\dot{x}) \qquad \Longrightarrow \qquad x(t) = \mathcal{X}(t,x_0,\dot{x}_0) \; . \tag{2}$ From this form it is apparent that solutions depend upon two pieces of initial value data: $x_0 = x(0)$ and $\dot{x}_0 = \dot{x}(0)$.

The fact that solutions require two pieces of initial value data means that there must be two canonical coordinates, $X$ and $P$. They are traditionally taken to be, $X \equiv x \qquad {\rm and} \qquad P \equiv \frac{\partial L}{\partial \dot{x}} \; . \tag{3}$ The assumption of nondegeneracy implies one can invert the phase space transformation (3) to solve for $\dot{x}$ in terms of $X$ and $P$. That is, there exists a velocity $V(X,P)$ such that, $\frac{\partial L}{\partial \dot{x}} \Biggl\vert_{x = X \atop \dot{x} = V} = P \; . \tag{4}$

The canonical Hamiltonian is obtained by Legendre transforming on $\dot{x}$, \begin{align} H(X,P) & \equiv P \dot{x} - L \; , \\ & = P V(X,P) - L\Bigl(X,V(X,P)\Bigr) \; . \end{align} It is easy to check that the canonical evolution equations reproduce the inverse phase space transformation (4) and the Euler-Lagrange equation (1), \begin{align} \dot{X} & \equiv \frac{\partial H}{\partial P} = V + P \frac{\partial V}{ \partial P} - \frac{\partial L}{\partial \dot{x}} \frac{\partial V}{\partial P} = V \; , \\ \dot{P} & \equiv -\frac{\partial H}{\partial X} = -P \frac{\partial V}{ \partial X} + \frac{\partial L}{\partial x} + \frac{\partial L}{\partial \dot{x}} \frac{\partial V}{\partial X} = \frac{\partial L}{\partial x} \; . \end{align} This is the meaning of the statement, the Hamiltonian generates time evolution. When the Lagrangian has no explicit time dependence, $H$ is also the associated conserved quantity. Hence it possesses the key properties physicists want for the energy, and is unique up to canonical transformations.

A familiar example is the simple harmonic oscillator of mass $m$ and frequency $\omega$ whose Lagrangian is, $L = \frac12 m \dot{x}^2 - \frac12 m \omega^2 x^2 \; . \tag{5}$ The equation of motion and its general solution are, $\ddot{x}(t) = -\omega^2 x(t) \qquad \Longrightarrow \qquad x(t) = x_0 \cos(\omega t) + \frac{ \dot{x}_0}{\omega} \sin(\omega t) \; . \tag{6}$ The canonical variables for this system are, $X = x \qquad {\rm and} \qquad P = m \dot{x} \qquad \Longrightarrow \qquad V(X,P) = \frac{P}{m} \; .$ And the Hamiltonian is, $H = \frac{P^2}{2 m} + \frac12 m \omega^2 X^2 \; .$ Because it is quadratic in both $X$ and $P$, the Hamiltonian $H(X,P)$ is bounded below by zero.

Now consider a system whose Lagrangian $L(x,\dot{x},\ddot{x})$ depends nondegenerately upon $\ddot{x}$. The Euler-Lagrange equation is, $\frac{\partial L}{\partial x} - \frac{d}{dt} \frac{\partial L}{\partial \dot{x}} + \frac{d^2}{dt^2} \frac{\partial L}{\partial \ddot{x}} = 0 \; . \tag{7}$ Now nondegeneracy means $\frac{\partial^2 L}{\partial \ddot{x}^2} \neq 0$, which implies that the Euler-Lagrange equation (7) can be cast in a form radically different from Newton's, $\ddddot{x} = \mathcal{F}(x,\dot{x},\ddot{x},\dddot{x}) \qquad \Longrightarrow \qquad x(t) = \mathcal{X}(t,x_0,\dot{x}_0,\ddot{x}_0, \dddot{x}_0) \; .$

Because solutions now depend upon four pieces of initial value data there must be four canonical coordinates. Ostrogradsky's choices for these are, \begin{align} X_1 \equiv x \qquad & , \qquad P_1 \equiv \frac{\partial L}{\partial \dot{x}} - \frac{d}{dt} \frac{\partial L}{\partial \ddot{x}} \; , \tag{8} \\ X_2 \equiv \dot{x} \qquad & , \qquad P_2 \equiv \frac{\partial L}{\partial \ddot{x}} \; . \tag{9} \end{align} The assumption of nondegeneracy implies one can invert the phase space transformation (8)-(9) to solve for $\ddot{x}$ in terms of $X_1$, $X_2$ and $P_2$. That is, there exists an acceleration $A(X_1,X_2,P_2)$ such that, $\frac{\partial L}{\partial \ddot{x}} \Biggl\vert_{{x = X_1 \atop \dot{x} = X_2} \atop \ddot{x} = A} = P_2 \; . \tag{10}$ Note that the acceleration $A(X_1,X_2,P_2)$ does nor depend upon $P_1$. The momentum $P_1$ is only needed for the third time derivative.

Ostrogradsky's Hamiltonian is obtained by Legendre transforming on $\dot{x} = x^{(1)}$ and $\ddot{x} = x^{(2)}$, \begin{align} H(X_1,X_2,P_1,P_2) &\equiv \sum_{i=1}^2 P_i x^{(i)} - L \; , \\ & = P_1 X_2 + P_2 A(X_1,X_2,P_2) - L\Bigl(X_1,X_2,A(X_1,X_2,P_2)\Bigr) \; . \tag{11} \end{align} The time evolution equations are those suggested by the notation, $\dot{X_i} \equiv \frac{\partial H}{\partial P_i} \qquad {\rm and} \qquad \dot{P}_i \equiv - \frac{\partial H}{\partial X_i} \; .$ To check that they generate time evolution, note first that the evolution equation for $X_1$ is, $\dot{X}_1 = \frac{\partial H}{\partial P_1} = X_2 \; .$ Of course this reproduces the phase space transformation $\dot{x} = X_2$ in (9). The evolution equation for $X_2$ similarly reproduces (10), $\dot{X}_2 = \frac{\partial H}{\partial P_2} = A + P_2 \frac{\partial A}{ \partial P_2} - \frac{\partial L}{\partial \ddot{x}} \frac{\partial A}{\partial P_2} = A \; .$ The phase space transformation $P_1 = \frac{\partial L}{\partial\dot{x}} - \frac{d}{dt} \frac{\partial L}{\partial \ddot{x}}$ (8) comes from the evolution equation for $P_2$, $\dot{P}_2 = -\frac{\partial H}{\partial X_2} = -P_1 - P_2 \frac{\partial A}{ \partial X_2} + \frac{\partial L}{\partial \dot{x}} + \frac{\partial L}{ \partial \ddot{x}} \frac{\partial A}{\partial X_2} = -P_1 + \frac{\partial L}{ \partial \dot{x}} \; .$ And the Euler-Lagrange equation (7) follows from the evolution equation for $P_1$, $\dot{P}_1 = -\frac{\partial H}{\partial X_1} = -P_2 \frac{\partial A}{\partial X_1} + \frac{\partial L}{\partial x} + \frac{\partial L}{\partial \ddot{x}} \frac{\partial A}{\partial X_1} = \frac{\partial L}{\partial x} \; .$ Hence Ostrogradsky's Hamiltonian generates time evolution. When the Lagrangian contains no explicit dependence upon time it is also the conserved Noether current. It is therefore the energy, again up to canonical transformation.

Ostrogradsky's Hamiltonian (11) is linear in the canonical momentum $P_1$, which means that no system of this form can be stable. In fact, there is not even any barrier to decay. Note the power and generality of the result: it applies to every Lagrangian $L(x,\dot{x},\ddot{x})$ which depends nondegenerately upon $\ddot{x}$, independent of the details. The only assumption is nondegeneracy, and that simply means one cannot eliminate $\ddot{x}$ by partial integration.

It is useful to consider a higher derivative example which depends upon a dimensionless parameter $\epsilon$ that quantifies its deviation from the simple harmonic oscillator (5), $L = -\frac{\epsilon m}{2 \omega^2} \ddot{x}^2 + \frac{m}2 \dot{x}^2 - \frac{m\omega^2}2 x^2 \; . \tag{12}$ The Euler-Lagrange equation and its general solution are, \begin{align} 0 & = -m \Bigl[ \frac{\epsilon}{\omega^2} \ddddot{x} + \ddot{x} + \omega^2 x\Bigr] \; , \tag{13} \\ x(t) & = C_+ \cos(k_+ t) + S_+ \sin(k_+ t) + C_- \cos(k_- t) + S_- \sin(k_- t) \; . \tag{14} \qquad \end{align} Here the two frequencies are, $k_{\pm} \equiv \omega \sqrt{ \frac{1 \mp \sqrt{1 \!-\! 4 \epsilon}}{2 \epsilon} } \; , \tag{15}$ and the constants $C_{\pm}$ and $S_{\pm}$ are functions of the initial value data, \begin{align} C_+ = \frac{k_-^2 x_0 \!+\! \ddot{x}_0}{k_-^2 \!-\! k_+^2} \qquad & , \qquad S_+ = \frac{k_-^2 \dot{x}_0 \!+\! \dddot{x}_0}{k_+ (k_-^2 \!-\! k_+^2)} \; , \\ C_- = \frac{k_+^2 x_0 \!+\! \ddot{x}_0}{k_+^2 \!-\! k_-^2} \qquad & , \qquad S_- = \frac{k_+^2 \dot{x}_0 \!+\! \dddot{x}_0}{k_- (k_+^2 \!-\! k_-^2)} \; . \end{align} For this model Ostrogradsky's two conjugate momenta (8)-(9) are, \begin{align} P_1 = m \dot{x} + \frac{\epsilon m}{\omega^2} \dddot{x} \qquad & \Leftrightarrow \qquad \dddot{x} = \frac{\omega^2 P_1 \!-\! m \omega^2 X_2}{\epsilon m} \; , \\ P_2 = - \frac{\epsilon m}{\omega^2} \ddot{x} \qquad & \Leftrightarrow \qquad \ddot{x} \equiv A = -\frac{\omega^2 P_2}{\epsilon m} \; . \end{align} The Hamiltonian can be expressed alternatively in terms of canonical variables, configuration space variables, or the constants $C_{\pm}$ and $S_{\pm}$, \begin{align} H & = P_1 X_2 - \frac{\omega^2}{2 \epsilon m} P_2^2 - \frac{m}2 X_2^2 + \frac{m\omega^2}2 X_1^2 \; , \tag{16} \\ & = \frac{\epsilon m}{\omega^2} \dot{x} \dddot{x} - \frac{\epsilon m}{2 \omega^2} \ddot{x}^2 + \frac{m}2 \dot{x}^2 + \frac{m \omega^2}2 x^2 \; , \tag{17} \\ & = \frac{m}2 \sqrt{1 \!-\! 4 \epsilon} \, k_+^2 (C_+^2 \!+\! S_+^2) - \frac{m}2 \sqrt{1 \!-\! 4 \epsilon} \, k_-^2 (C_-^2 \!+\! S_-^2) \; . \tag{18} \end{align} The last form (18) makes it clear that the $+$ modes carry positive energy and the $-$ modes carry negative energy.

Ostrogradsky's Construction for $N$ Derivatives

Consider a Lagrangian $L\left(x,\dot{x},\dots,x^{(N)}\right)$ which depends upon the first $N$ derivatives of $x(t)$. If this Lagrangian depends nondegenerately upon the $N$-th derivative $x^{(N)}$ then the Euler-Lagrange equation is linear in the $2N$-th derivative $x^{(2N)}$, $\sum_{i=0}^N \left(-{d \over dt}\right)^i {\partial L \over \partial x^{(i)}} = 0 \; . \tag{19}$ The canonical phase space must therefore possess $2N$ coordinates which Ostrogradsky chooses to be, $X_i \equiv x^{(i-1)} \qquad {\rm and} \qquad P_i \equiv \sum_{j=i}^N \Bigl(- \frac{d}{dt}\Bigr)^{j-i} \frac{\partial L}{\partial x^{(j)}} \; . \tag{20}$ Nondegeneracy means one can solve for $x^{(N)}$ in terms of $P_N$ and the $X_i$'s. That is, there exists a function $\mathcal{A}(X_1,\ldots,X_N,P_N)$ such that, $\frac{\partial L}{\partial q^{(N)}} \Biggl\vert_{x^{(i-1)} = X_i \atop x^{(N)} = \mathcal{A}} = P_N \; . \tag{21}$

For general $N$ Ostrogradsky's Hamiltonian takes the form, \begin{align} H & \equiv \sum_{i=1}^N P_i x^{(i)} - L \; , \\ & = P_1 X_2 + P_2 X_3 + \cdots + P_{N-1} X_N + P_N \mathcal{A} - L\Bigl(X_1,\ldots,X_N,\mathcal{A}\Bigr) \; . \tag{22} \end{align} The evolution equations are, $\dot{X}_i \equiv \frac{\partial H}{\partial P_i} \qquad {\rm and} \qquad \dot{P}_i \equiv -\frac{\partial H}{\partial X_i} \; .$ It is simple to check that these evolution equations reproduce the canonical relations (20) and the Euler-Lagrange equation (19). The first $(N-1)$ equations for $X_i$ verify the definition of $X_{i+1}$, $i=1,\dots,(N-1) \qquad \Longrightarrow \qquad \dot{X}_i = X_{i+1} \; .$ The evolution equation for $X_N$ is similar, $\dot{X}_N = \mathcal{A} + P_N \frac{\partial \mathcal{A}}{\partial P_N} - \frac{\partial L}{\partial x^{(N)}} \frac{\partial \mathcal{A}}{ \partial P_N} = \mathcal{A} \; .$ The last $(N-1)$ equations for $P_i$ reproduce the definition of $P_{i-1}$, \begin{align} i=2,\dots,N \quad \Longrightarrow \quad \dot{P}_i & = -P_{i-1} - P_N \frac{\partial \mathcal{A}}{\partial X_i} + \frac{\partial L}{\partial x^{(i-1)}} + \frac{\partial L}{\partial x^{(N)}} \frac{\partial \mathcal{A}}{\partial X_i} \; , \qquad \\ & = -P_{i-1} + \frac{\partial L}{\partial x^{(i-1)}} \; . \end{align} And the evolution equation for $P_1$ gives the Euler-Lagrange equation (19), $\dot{P}_1 = -P_N \frac{\partial \mathcal{A}}{\partial X_1} + \frac{\partial L}{\partial x} + \frac{\partial L}{\partial x^{(N)}} \frac{\partial \mathcal{A}}{\partial X_1} = \frac{\partial L}{\partial x} \; .$ Hence (22) generates time evolution. It is also the Noether current for the case where the Lagrangian contains no explicit time dependence. The Hamiltonian (22) is therefore what any physicist would call the energy, up to canonical transformation.

The Hamiltonian (22) is linear in $P_1, P_2, \ldots P_{N-1}$. Only with respect to $P_N$ might it be bounded from below. For large $N$ the fraction of linear directions approaches $\frac12$, so adding more higher derivatives makes the instability worse rather than better.

Nature of the Instability

Ostrogradsky's result implies that the Hamiltonian of a nondegenerate higher derivative theory is unbounded below, and also above. This section discusses the manner in which the instability manifests, and what it implies for fundamental theory. Six short subsections make the points:

1. The Ostrogradskian instability drives the dynamical variable to a special kind of time dependence, not a special numerical value.
2. The same Ostrogradskian dynamical variable carries both positive and negative energy creation and annihilation operators.
3. If a system which suffers from the Ostrogradskian instability interacts, then the empty state can decay into a collection of positive and negative energy excitations.
4. If a system which suffers from the Ostrogradskian instability is a continuum field theory, the vast entropy at infinite 3-momentum will make the decay instantaneous.
5. For interacting systems which suffer from the Ostrogradskian instability, degrees of freedom with large 3-momentum do not decouple from low energy physics.
6. The imposition of a single, global constraint on the energy functional does not ameliorate the Ostrogradskian instability.

Kinetic Instability

Physicists are familiar with instabilities of the potential energy. In this case energy is released as the dynamical variable approaches some special value. The Ostrogradskian instability is instead a problem with the kinetic energy, and it manifests by the dynamical variable developing a special time dependence. Checking that the energy is bounded below for constant values of the dynamical variable in no way establishes that a system is free of the Ostrogradskian instability. Consider, for example, the higher derivative oscillator (12). Expression (17) shows that its energy is bounded below by zero for any constant value of $x(t)$. Negative energies are attained by making $\ddot{x}(t)$ large and/or making $\dddot{x}(t)$ large while keeping the combination $\dot{x}(t) \!+\! \frac{\epsilon}{\omega^2}\dddot{x}(t)$ fixed.

Double Duty for Dynamical Variables

Physicists are used to resolving linearized dynamical variables into creation and annihilation operators. For the harmonic oscillator solution (6) this is done by using the Euler relation to identify a lowering operator proportional to $e^{-i\omega t}$ and a raising operator proportional to $e^{i \omega t}$, $x_0 \cos(\omega t) + \frac{ \dot{x}_0}{\omega} \sin(\omega t) = \frac12 \Bigl[ x_0 + \frac{i}{\omega} \dot{x}_0\Bigr] e^{-i \omega t} + \frac12 \Bigl[ x_0 - \frac{i}{\omega} \dot{x}_0\Bigr] e^{i \omega t} \; .$ The usual rule is that each dynamical variable harbors either zero or one set of creation and annihilation operators at linearized order. From expression (14) one can see that the same higher derivative dynamical variable carries both positive and negative energy creation and annihilation operators. This means that local interactions which involve the dynamical variable necessarily couple the two sectors.

The Vacuum Can Decay

Now consider an interacting, continuum field theory which possesses the Ostrogradskian instability. In particular consider its likely particle spectrum about some “empty” solution in which the field is constant. Because the Hamiltonian is linear in all but one of the conjugate momenta it is possible to arbitrarily increase or decrease the energy by moving different directions in phase space. Hence there must be both positive energy and negative energy particles -- just as there are in the higher derivative oscillator (12). As in that point particle model, the same continuum field must carry the creation and annihilation operators of both the positive and the negative energy particles. If the theory is interacting at all -- that is, if its Lagrangian contains a higher than quadratic power of the field -- then there will be interactions between positive and negative energy particles. Depending upon the interaction, the empty state can decay into some collection of positive and negative energy particles.

Entropy Drives Vacuum Decay

Recall the reason that excited states of atoms decay in nature. It is certainly not to reduce the energy of the full system -- including the interaction with electromagnetism -- but rather to redistribute the constant total energy into the largest possible class of states. There is one way for the atom not to decay, compared with an infinite number of ways the atom can decay and emit one or more recoil photons. Note also that explicit computations of the decay time employ vacuum fluctuations of the electromagnetic field to provide the necessary perturbation.

Atomic decays have just the fixed energy difference between the two states to apportion, so they are chiefly driven by the arbitrary directions which can be taken by the decay products. In contrast, the decay of an interacting, nondegenerate higher derivative field theory can involve particles of {\it any} energy, as long as the total sums to zero. So one should think of the decay rate as having the same sort of angular factors as an atomic decay at some fixed energy, followed by one or more integrals -- all the way to infinity --- over the magnitudes of the various energies. The volume of phase space is so large that these integrations cause the decay to be instantaneous. Indeed, the only way people derive finite decay rates for particles with a kinetic instability is by cutting off the phase space at some point, in which case the rate is dominated by the cutoff, for example see Cline et al., 2004. Such a cutoff might make sense if the kinetic instability appeared in some nonlocal effective field theory, but it has no place in fundamental physics.

Note that the decay does not just happen once. It is even more entropicly favored for there to be two decays, and better yet for more. In fact the system instantly evaporates into a maelstrom of positive and negative energy particles. Whether or not such a state has a proper mathematical representation, it certainly does not describe the universe of human experience in which all particles have positive energy and empty space remains empty.

Note also that this conclusion only follows if the higher derivative theory possesses both interactions and continuum particles. The point particle oscillator (12) has no interactions, so its negative energy degree of freedom is unobservable. However, it is conceivable that this higher derivative oscillator could be coupled to a discrete system without engendering any instability. The feature which drives explosive vacuum decay is the vast entropy of phase space. Without that it becomes an open question whether or not there is anything wrong with a higher derivative theory. Of course the physical universe seems to be described by continuum field theory down to at least $2.8 \times 10^{-19}$ meters (Bourlikov, 2000), and any observable degree of freedom must interact, or else it could not be observed, so these seem to be safe assumptions.

Large $\Vert \vec{k}\Vert$ Modes Do Not Decouple

Physicists are used to ignoring very high energy modes, except for renormalizations of low energy parameters. This procedure is quite correct for positive energy modes in a stable theory because exciting a mode requires energy which must be drawn from de-exciting other modes, and any given state only has some fixed amount of energy. However, that justification fails for a theory which suffers from the Ostrogradskian instability because even a very high (positive or negative) energy mode can be excited by also exciting modes with the opposite energy. Instead of these large $k$ modes decoupling, they couple ever more strongly as $k$ grows, because more and more ways open up to balance its energy by exciting lower modes of the opposite sign.

Constraints on $H$ Accomplish Nothing

It is sometimes imagined that the energy of a higher derivative theory decays with time. That is not true. Provided one is dealing with a complete system, and provided there is no external time dependence, the energy of a higher derivative system is conserved, just as it would be under those conditions for a lower derivative theory. This conservation is apparent for the higher derivative oscillator (12) from expression (18).

The physical problem with nondegenerate higher derivative theories is not that their energies decay to lower and lower values. The problem is rather that certain sectors of the theory become arbitrarily highly excited when one is dealing with an interacting, continuum field theory which has nondegenerate higher derivatives. For example, Boulware et al., 2000 showed that the energy is zero for any asymptotically flat solution of the higher derivative field equations derived from the Lagrangian, $\mathcal{L} = -\frac14 \alpha C_{\mu\nu\rho\sigma} C^{\mu\nu\rho\sigma} \sqrt{-g} -\frac14 \beta R^2 \sqrt{-g} \; , \tag{23}$ where $C_{\mu\nu\rho\sigma}$ is the Weyl tensor and $R$ is the Ricci scalar. However, this model is still unstable for $\alpha \neq 0$, as its creators realized. (It is also worth noting that the requirement of asymptotic flatness in this model would preclude the response to normal matter, and that imposing the correct asymptotic condition gives rise to nonzero energy Deser and Tekin, 2002).

Quantization

Quantization is very important to understanding the Ostrogradskian instability because 0-point fluctuations provide the perturbations needed to ensure that the potential for vacuum decay is actually realized. However, a quantum higher derivative system has some peculiarities. For example, it is obvious from relations (8)-(9) that position and velocity commute! Further, the wave function of a higher derivative theory depends upon position and velocity. This section argues first that the classical instability survives canonical quantization. After presenting a worked-out example, the curious noncanonical quantizations which sometimes appear in the literature are discussed.

A Large Phase Space Instability

It is often imagined that quantization might protect a higher derivative system against the Ostrogradskian instability the same way that quantization prevents the collapse of atoms coupled to electromagnetism. This is a failure to understand correspondence limits. In the Heisenberg picture the equations of classical mechanics are identical to those of quantum mechanics. It also means the very same thing to solve these equations: one expresses the dynamical variable in terms of time and the allowed initial value data, as in expressions (6) and (14). The only difference between classical and quantum mechanics is that the classical initial value data are numbers which can take any value whereas the quantum initial value data include noncommuting conjugate operators which obey the Uncertainty Principle. The only classical phenomena that can be affected by quantization are those whose realization requires localizing conjugate variables to some volume of the classical phase space smaller than $\hbar$. So quantum atoms are stable because localizing the electrons too near the nucleus necessarily induces a large kinetic energy.

In contrast, the Ostrogradskian instability derives from the fact that $P_1 X_2$ can be made arbitrarily negative by taking $P_1$ either very negative, for positive $X_2$, or else very positive, for negative $X_2$. This covers essentially half the classical phase space! Further, the variables $X_2$ and $P_1$ commute with one another in Ostrogradskian quantum mechanics. So there is no reason to expect that the Ostrogradskian instability is unaffected by quantization.

Quantum Higher Derivative Oscillator

Consider the second derivative oscillator (12) discussed in the section on “Ostrogradsky’s construction for two derivatives”. There can be no ground state in the presence of the Ostrogradskian instability but one might define an “empty” state wavefunction, $\Omega(X_1,X_2)$ which has the minimum excitation in both the positive and negative energy degrees of freedom. The procedure for doing this is simple: first identify the positive and negative energy lowering operators $\alpha_{\pm}$, and then solve the equations, $\alpha_+ \vert \Omega \rangle = 0 = \alpha_- \vert \Omega \rangle \; . \tag{24}$ One can recognize the raising and lowering operators by expressing the general solution (14) in terms of exponentials, $x(t) = \frac12 (C_+ \!+\! i S_+) e^{-ik_+ t} + \frac12 (C_+ \!-\! i S_+) e^{ik_+ t} + \frac12 (C_- \!+\! i S_-) e^{-ik_-t } + \frac12 (C_- \!-\! i S_-) e^{ik_- t} \; .$ Recall that the $k_+$ mode carries positive energy, so its lowering operator must be proportional to the $e^{-ik_+ t}$ term, \begin{align} \alpha_+ & \sim C_+ + i S_+ \; , \\ & \sim \frac{m k_+}2 \Bigl(1 \!+\! \sqrt{1 \!-\! 4\epsilon}\Bigr) X_1 + i P_1 - k_+ P_2 - \frac{i m}2 \Bigl(1 \!-\! \sqrt{1 \!-\! 4\epsilon}\Bigr) X_2 \; . \end{align} The $k_-$ mode carries negative energy, so its lowering operator must be proportional to the $e^{+i k_- t}$ term, \begin{align} \alpha_- & \sim C_- - i S_- \; , \\ & \sim \frac{m k_-}2 \Bigl(1 \!-\! \sqrt{1 \!-\! 4\epsilon}\Bigr) X_1 - i P_1 - k_- P_2 + \frac{i m}2 \Bigl(1 \!+\! \sqrt{1 \!-\! 4\epsilon}\Bigr) X_2 \; . \end{align} Writing $P_i = -i\hbar \frac{\partial}{\partial X_i}$ reveals that the unique solution to (24) has the form, $\Omega(X_1,X_2) = N \exp\Biggl[-\frac{m \sqrt{1 \!-\! 4\epsilon}}{2 \hbar(k_+ \!+\! k_-)} \Bigl(k_+ k_- X_1^2 + X_2^2\Bigr) - \frac{i \sqrt{\epsilon} m}{\hbar} X_1 X_2\Biggr] \; . \tag{25}$

The empty wave function (25) is obviously normalizable, so it gives a state of the quantum system. One can build a complete set of normalized stationary states by acting arbitrary numbers of $+$ and $-$ raising operators on it, $\vert N_+ , N_-\rangle \equiv \frac{(\alpha_+^{\dagger})^N_+}{\sqrt{N_+ !}} \frac{(\alpha_-^{\dagger})^N_-}{\sqrt{N_- !}} \vert \Omega \rangle \; .$ On this space of states the Hamiltonian operator is unbounded below, just as in the classical theory, $H \vert N_+ , N_- \rangle = \hbar \Bigl(N_+ k_+ - N_- k_-\Bigr) \vert N_+ , N_- \rangle \; .$ This is the correct way to quantize a higher derivative theory. One evidence of this fact is that classical configurations of negative energy correspond to quantum negative energy states.

Unitarity versus Instability

Particle physicists who quantize higher derivative theories do not typically recognize a problem with stability; they instead discuss a breakdown of unitarity, for example see Stelle, 1977. This is accomplished by regarding the negative energy lowering operator as a positive energy raising operator. So one defines a “ground state” $\vert \overline{\Omega} \rangle$ which obeys the equations, $\alpha_+ \vert \overline{\Omega} \rangle = 0 = \alpha_-^{\dagger} \vert \overline{\Omega} \rangle \; . \tag{26}$ The unique wave function which solves these equations is, $\overline{\Omega}(X_1,X_2) = N \exp\Biggl[-\frac{m \sqrt{1 \!-\! 4\epsilon}}{2 \hbar (k_- \!-\! k_+)} \Bigl(k_+ k_- X_1^2 - X_2^2\Bigr) + \frac{i \sqrt{\epsilon} m}{\hbar} X_1 X_2 \Biggr] \; . \tag{27}$ The wave function (27) is not normalizable, so it does not correspond to a state of the quantum system Hawking and Hertog, 2002. However, particle physicists define a formal “space of states” based upon $\vert\overline{\Omega} \rangle$, $\vert \overline{N_+ , N_-}\rangle \equiv \frac{(\alpha_+^{\dagger})^{N_+}}{ \sqrt{N_+ !}} \frac{(\alpha_-)^{N_-}}{\sqrt{N_- !}} \vert \overline{\Omega} \rangle \; .$ Although these wave functions are no more normalizable than $\overline{\Omega}(X_1,X_2)$, they are all positive energy eigenfunctions, $H \vert \overline{N_+ , N_-} \rangle = \hbar \Bigl(N_+ k_+ + N_- k_-\Bigr) \vert \overline{N_+ , N_-} \rangle \; .$

The problem with unitarity emerges because $\vert \overline{\Omega}\rangle$ is defined to have unit norm, but the commutation relations are unchanged, $[\alpha_+,\alpha_+^{\dagger}] = 1 = [\alpha_-,\alpha_-^{\dagger}] \; .$ Hence the norm of any state with odd $N_-$ is negative. The first of these negative norm states is, $\langle \overline{0,1} \vert \overline{0,1}\rangle = \langle \overline{\Omega} \vert \alpha_-^{\dagger} \alpha_- \vert \overline{ \Omega} \rangle = - \langle \overline{\Omega} \vert \overline{\Omega} \rangle \; .$ The next step is to invoke the probabilistic interpretation of quantum mechanics which requires norms to be positive because probabilities are. Therefore, the negative norm states must be excised from the space of states. However, doing that results in a nonunitary S-matrix because scattering processes inevitably mix positive and negative norm states, just as the correctly-quantized, indefinite-energy theory allows processes which mix positive and negative energy particles.

It is important to note that the potential for invoking noncanonical quantization schemes to change the range of allowed energies is present even in the usual, first derivative systems. The Schrodinger equation $H \psi(X) = E \psi(X)$ is a second order differential equation, which possesses two linearly independent solutions for every value of the energy $E$. It is only by insisting upon normalizable wave functions that quantized energies emerge. Many other peculiar things happen if one abandons normalizability (Woodard, 1993; Tsamis and Woodard, 1987). In particular, the Correspondence Principle fails, so that taking $\hbar$ to zero gives a different classical system from the one which originally motivated the analysis. That is the case for $PT$-symmetric quantizations of higher derivative systems (Bender and Mannheim, 2008; Mostafazadeh, 2010).

Degeneracy

The only way anyone has ever found to avoid the Ostrogradskian instability is by violating the assumption of nondegeneracy upon which it is based. This section discusses three ways this can happen: through partial integration, through gauge invariance, and by imposing constraints by fiat Chen et al., 2013.

Trivial Degeneracy

The simplest form of degeneracy derives from adding a total derivative to a first order system. Examples include the Hilbert action of general relativity, Lovelock gravity (Lovelock, 1971) and Galileons (Nicolis et al., 2009; Rham and Tolley, 2010). In that case one simply performs a partial integration, and discards the surface term to obtain a Lagrangian which contains only first time derivatives. For example, the 3rd Lagrangian for a scalar Galileon $\pi(t,\vec{x})$ reduces to first order form as, \begin{align} \partial_{\mu} \pi \partial^{\mu} \pi \partial^2 \pi & = \frac{\partial}{\partial t} \Bigl[\frac13 \dot{\pi}^3 \!-\! \dot{\pi} \vec{\nabla} \pi \!\cdot\! \vec{\nabla} \pi\Bigr] + 2 \dot{\pi} \vec{\nabla} \pi \!\cdot\! \vec{\nabla} \dot{\pi} + \nabla^2 \pi \vec{\nabla} \pi \!\cdot\! \vec{\nabla} \pi \; , \qquad \\ & \longrightarrow 2 \dot{\pi} \vec{\nabla} \pi \!\cdot\! \vec{\nabla} \dot{\pi} + \nabla^2 \pi \vec{\nabla} \pi \!\cdot\! \vec{\nabla} \pi \; . \end{align} Note that it is only necessary to eliminate higher {\it time} derivatives; there is no problem if the Lagrangian contains higher spatial derivatives, or mixed first time and space derivatives.

Gauge Degeneracy

All theories which possess continuous symmetries are degenerate, irrespective of whether or not they possess higher derivatives. A familiar example is the relativistic point particle, whose dynamical variable is $X^{\mu}(\tau)$ and whose Lagrangian is, $L = -m \sqrt{-\eta_{\mu\nu} \dot{X}^{\mu} \dot{X}^{\nu}} \; .$ The conjugate momentum is, $P_{\mu} \equiv \frac{m \dot{X}_{\mu}}{\sqrt{-\dot{X}^2}} \; . \tag{28}$ One cannot solve (28) for $\dot{X}^{\mu}$ in terms of $X^{\mu}$ and $P_{\mu}$ because the equation is homogeneous of degree zero. The continuous symmetry associated with this degeneracy is invariance under changes of the parameter $\tau \longrightarrow \tau'$, $X^{\mu}(\tau) \longrightarrow X^{\prime \mu}(\tau) \equiv X^{\mu}\Bigl({\tau' }^{-1}(\tau)\Bigr) \; .$

The cure for symmetry-induced degeneracy is simply to fix the symmetry by imposing gauge conditions. Then the gauge-fixed Lagrangian should no longer be degenerate in terms of the remaining variables. For example, one might fix the parameter $\tau$ to obey $\tau = X^0(\tau)$. In that case the gauge-fixed particle Lagrangian is, $L_{\rm GF} = -m \sqrt{1 - \dot{\vec{X}} \cdot \dot{\vec{X}} } \; ,$ and the relations for the momenta are simple to invert, $P_i \equiv \frac{m \dot{X}_i}{\sqrt{1 - \dot{\vec{X}} \cdot \dot{\vec{X}} }} \qquad \Longleftrightarrow \qquad \dot{X}^i = \frac{P^i}{\sqrt{m^2 + \vec{P} \cdot \vec{P}}} \; .$

When a continuous symmetry is used to eliminate a dynamical variable, the equation of motion of this variable typically becomes a constraint. For symmetries enforced by means of a compensating field -- such as making the Hilbert action local Lorentz invariant using the antisymmetric components of the vierbein (Woodard, 1984), or Weyl invariant using a scalar (Tsamis and Woodard, 1986) -- the associated constraints are tautologies of the form $0 = 0$. Sometimes the constraints are nontrivial, but implied by the equations of motion. An example of this kind is the relativistic particle considered above. In synchronous gauge ($\tau = X^0(\tau)$) the equation of the gauge-fixed zero-component implies that the Hamiltonian is conserved, $\frac{d}{d\tau} \Biggl( \frac{m \dot{X}_0}{\sqrt{-\eta_{\mu\nu} \dot{X}^{\mu} \dot{X}^{\nu}}} \Biggr) = 0 \longrightarrow \frac{d}{dt} \Bigl( \sqrt{m^2 + \vec{P} \!\cdot\! \vec{P} } \Bigr) = 0 \; .$ And sometimes the constraints give nontrivial relations between the canonical variables that generate residual, time-independent symmetries. In this case another degree of freedom can be removed. An example of this kind of constraint is Gauss' Law in temporal gauge electrodynamics.

When constraints of the third type are present one must check whether or not they affect the instability. This obviously depends on the particular model being studied but a necessary condition for avoiding the Ostrogradskian instability is that the number of gauge constraints must equal or exceed the number of unstable directions in the canonical phase space. Because the number of constraints for any given symmetry is fixed, whereas the number of unstable directions increases with the number of higher derivatives, it follows that gauge constraints can at best avoid instability for some fixed number of higher derivatives.

A good example of gauge degeneracy is provided by the quadratic curvature model (23) which was exhibited at the end of the section on “Nature of the instability” to show the irrelevance of a global constraint on the Hamiltonian. As long as $\alpha$ and $\beta$ are both nonzero, there are six independent, higher derivative momenta at each space point, whereas there are only four local constraints. If $\beta = 0$ the model acquires a new local symmetry -- Weyl invariance -- which adds another local constraint. Hence there are either one or two unconstrained instabilities per space point for $\alpha \neq 0$. There are an infinite number of space points, so the addition of a single, global constraint does not change anything.

The case of $\alpha = 0$ is special. If $\beta$ has the right sign the resulting model has long been known to have positive energy (Starobinsky, 1980; Strominger, 1984). This result in no way contradicts the previous analysis. When $\alpha = 0$, the terms which carry second derivatives are contracted in such way that only a single component of the metric carries higher derivatives. So the counting is one unstable direction per space point versus four local constraints, which means the constraints can prevent the Ostrogradskian instability.

Imposed Degeneracy

Many attempts to evade the Ostrogradskian instability are based on segregating higher derivatives to interaction terms so that the free theory possesses no extra solutions. This renders the instability invisible to perturbative scrutiny but does not avoid it. One can see from the construction of Ostrogradsky that the sole assumption needed to derive the instability is nondegeneracy, irrespective of how one organizes any approximation technique. On the other hand, there is a way of imposing constraints so as to make the theory agree with its perturbative development. When this is done there are no more higher derivative degrees of freedom, but this constrained version of the theory cannot serve to define an acceptable model unless the perturbative solution converges.

The technique is to regard higher derivative parts of the Euler-Lagrange equation as a perturbation and then use the unperturbed equation to reduce the order (Jaén et al., 1986). Of course this produces a remainder with even more higher derivatives, but this remainder is also higher order in perturbation theory. By iterating the procedure infinitely, and then neglecting the remainder, one obtains a lower order equation.

The technique can be illustrated for the higher derivative oscillator (12) by regarding the parameter $\epsilon$ as a coupling constant so that the Euler-Lagrange equation (13) takes the form, $\ddot{x} + \omega^2 x = -\epsilon \Bigl( \frac{d}{d \omega t}\Bigr)^2 \ddot{x} \equiv -\epsilon D^2 \ddot{x} \; . \tag{29}$ The first iteration gives, $\ddot{x} + \omega^2 x = + \epsilon \ddot{x} + \epsilon^2 D^4 \ddot{x} = -\epsilon \omega^2 x - \epsilon^2 (1 \!-\! D^2) D^2 \ddot{x} \; . \tag{30}$ After another iteration one obtains, \begin{align} \ddot{x} + \omega^2 x &= -\epsilon \Bigl[1 \!+\! \epsilon (1 \!+\! \epsilon)^2 (2 \!+\! \epsilon)\Bigr] \omega^2 x - \epsilon^4 \Bigl[(2 \!+\! \epsilon)(1\!+\! \epsilon) - (2\!+\!\epsilon) D^2 + D^4\Bigr] (1 \!-\! D^2) D^2 \ddot{x} \; . \qquad \tag{31} \end{align} Continuing in this fashion, and ignoring the remainder, gives, $\ddot{x} + k_+^2 x = 0 \; . \tag{32}$ From the full theory, the perturbative development has retained only the solution whose frequency is well behaved for $\epsilon \rightarrow 0$, $k_+^2 = \omega^2 \Bigl[1 + \epsilon + 2 \epsilon^2 + O(\epsilon^3) \Bigr] \; . \tag{33}$ It has discarded the solution whose frequency blows up as $\epsilon \rightarrow 0$, $k_-^2 = \omega^2 \Bigl[\frac1{\epsilon} - 1 - \epsilon - 2 \epsilon^2 + O(\epsilon^3)\Bigr] \; . \tag{34}$

The perturbative development (32) is what results if one changes the original theory by imposing the constraints, \begin{align} \ddot{q}(t) = - k_+^2 q(t) \qquad & \Longleftrightarrow \qquad P_2 = \frac{m}2 \Bigl(1 \!-\! \sqrt{1 \!-\! 4\epsilon}\Bigr) X_1 \; , \tag{35} \\ q^{(3)}(t) = - k_+^2 \dot{q}(t) \qquad & \Longleftrightarrow \qquad X_2 = \frac1{2 \epsilon m} \Bigl(1 \!-\! \sqrt{1 \!-\! 4\epsilon}\Bigr) P_1 \; . \tag{36} \end{align} Under these constraints the Hamiltonian becomes, $H_{\rm pert} = \sqrt{1 \!-\! 4\epsilon} \Bigl( \frac{m}2 X_2^2 + \frac{m k_+^2}2 X_1^2 \Bigr) \; ,$ which is that of a positive energy harmonic oscillator with mass $\sqrt{1 - 4\epsilon} m$ and frequency $k_+$. If the constraints (35)-(36) are imposed at one instant, they remain valid for all times as a consequence of the full equation of motion (13), so the constrained model is consistent. This is ultimately a consequence of the fact that, for this model, the perturbatve expansion converges. That is what ensures that the discarded remainder term really goes to zero when the expansion is carried to infinite order.

For nonlinear Euler-Lagrange equations it is more difficult to reach a second order form, but one can still do it. As before, the ultimate consistency of the reduced system depends upon the convergence of the perturbative expansion. For certain mechanical systems it does converge, for example, a higher derivative generalization of a particle moving in a uniform gravitational acceleration $g$ is, $L = \frac12 m \dot{x}^2 + m g x + \frac{\epsilon m}{6 g} x \ddot{x}^2 \qquad \Longrightarrow \qquad \ddot{x} = g + \frac{\epsilon}{6 g} \ddot{x}^2 + + \frac{\epsilon}{3 g} \frac{d^2}{dt^2} \Bigl( x \ddot{x} \Bigr) \; . \tag{37}$ Reducing to second order transforms the higher derivative corrections into a distortion of the acceleration, $\ddot{x} = \frac{g}{2 \epsilon} \Bigl[1 - \sqrt{1 \!-\! 2 \epsilon}\Bigr] = g \Bigl[1 + \frac12 \epsilon + \frac12 \epsilon^2 + O(\epsilon^3)\Bigr] \; . \tag{38}$ However, there are no known, interacting, $3+1$-dimensional field theories for which the perturbative expansion converges. Nor has anyone ever found a consistent way of imposing constraints which avoids the Ostrogradskian instability for an interacting, $(3 + 1)$-dimensional, higher derivative field theory.

Conclusions

Although it was not apparent in 1850, Ostrogradsky's theorem can today be recognized as the strongest restriction on what sorts of interacting local quantum field theories can describe fundamental physics. No symmetry principle has a broader scope or comparable power. Its applications include:

• Demonstrating that higher derivative counterterms cannot be a fundamental solution to the problem of quantum gravity (Woodard, 2009);
• Establishing $f(R)$ models as the only metric-based, local and potentially stable modifications of gravity (Woodard, 2007); and
• Discussing the problems of nonlocal models which can be viewed as the limits of an infinite sequence of higher derivatives (Eliezer and Woodard, 1989).

One should also note the recent generalization by Motohashi and Suyama of Ostrogradsky's result to Lagrangian-based systems (such as fermions) whose Euler-Lagrange equations involve an odd number of time derivatives Motohashi and Suyama, 2015.

Acknowledgements

It is a pleasure to acknowledge correspondence on this subject with S. Deser and B. P. Kosyakov. This work was partially supported by NSF grant PHY-1205591 and by the Institute for Fundamental Theory at the University of Florida.

References

• E. Straus, in Helle Zeit -- Dunkle Zeit (Europa Verlag A. G., Zurich, 1956) ed. C. Seelig, p. 72.
• M. Ostrogradsky, Mem. Ac. St. Petersbourg VI 4 (1850) 385.
• J. J. O'Connor and E. F. Robertson, “Mikhail Vasilevich Ostrogradski,” in MacTutor History of Mathematics, http://www-history.mcs.st-andrews.ac.uk/Biographies/Ostrogradski.html.
• J. M. Cline, S. Jeon and G. D. Moore, Phys. Rev. D 70 (2004) 043543, hep-ph/0311312.
• D. Bourilkov, Phys. Rev. D 62 (2000) 076005, hep-ph/0002172.
• D. G. Boulware, G .T. Horowitz and A. Strominger, Phys. Rev. Lett. 50 (1983) 1726.
• K. S. Stelle, Phys. Rev. D 16 (1977) 953.
• S. W. Hawking and T. Hertog, Phys. Rev. D 65 (2002) 103515, hep-th/0107088.
• R. P. Woodard, Class. Quant. Grav. 10 (1993) 483.
• N. C. Tsamis and R. P. Woodard, Phys. Rev. D 36 (1987) 3641.
• C. M. Bender and P. D. Mannheim, Phys. Rev. Lett. 100 (2008) 110402, arXiv:0706.0207.
• A. Mostafazadeh, Int. J. Geom. Meth. Mod. Phys. 7 (2010) 1191, arXiv:0810.5643.
• D. Lovelock, J. Math. Phys. 12 (1971) 498.
• A. Nicolis, R. Rattazzi and E. Trincherini, Phys. Rev. D 79 (2009) 064036, arXiv:0811.2197.
• C. de Rham and A. J. Tolley, JCAP 1005 (2010) 015, arXiv:1003.5917.
• R. P. Woodard, Phys. Lett. B 148 (1984) 440.
• N. C. Tsamis and R. P. Woodard, Ann. Phys. 168 (1986) 457.
• A. A. Starobinsky, Phys. Lett. B 91 (1980) 99.
• A. Strominger, Phys. Rev. D 30 (1984) 2257.
• X. Jaén, J. Llosa and A. Molina, Phys. Rev. D 34 (1986) 2302.
• R. P. Woodard, Rept. Prog. Phys. 72 (2009) 126002, arXiv:0907.4238.
• R. P. Woodard, Lect. Notes. Phys. 720 (2007) 403, astro-ph/0601672.
• D. A. Eliezer and R. P. Woodard, Nucl. Phys. B 325 (1989) 389.
• S. Deser and B. Tekin, Phys. Rev. Lett. 89 (2002), hep-th/0205318; S. Deser and B. Tekin, Phys. Rev. D 67 (2003), hep-th/0212292.
• T. J. Chen, M. Fasiello, E. A. Lim and A. J. Tolley, JCAP 1302 (2013) 042, arXiv:1209.0583.
• H. Motohashi and T. Suyama, Phys. Rev. D 91 (2015)085009, arXiv:1411.3721.