# The Inverted Pendulum as a Classical Analog of the EFT Paradigm

Martin Beneke  Matthias König  Martin Link  Physik Department T31, James-
Franck-Straße 1, Technische Universität München, D-85748 Garching, Germany

(August 26, 2023)

###### Abstract

The inverted pendulum is a mechanical system with a rapidly oscillating pivot
point. Using techniques similar in spirit to the methodology of effective
field theories, we derive an effective Lagrangian that allows for the
systematic computation of corrections to the so-called Kapitza equation. The
derivation of the effective potential of the system requires non-trivial
matching conditions, which need to be determined order by order in the power-
counting of the problem. The convergence behavior of the series is
investigated on the basis of high-order results obtained by this method.
TUM-HEP-1471/23

##  I The Inverted Pendulum

The Kapitza or inverted pendulum is a mechanical system consisting of an
inverted rigid pendulum with a rapidly oscillating pivot point P as depicted
in Fig. 1. The system is an interesting exercise for students of classical
mechanics due to a curious counter-intuitive phenomenon that occurs for a
certain range of parameters: If the frequency of the pivot oscillation is much
larger than the natural frequency of the unperturbed system and the amplitude
of said oscillation $h$ is smaller (within a certain range) than the length of
the pendulum $l$, the system exhibits a stable equilibrium in the upright
position Stephenson (1908); Kapitza (1951); Landau and Lifshitz (1960);
Butikov (2001).

Figure 1: Depiction of the inverted pendulum with shaft of length $l$ and
point mass $m$. The point O denotes the coordinate origin $y=0$ around which
the pivot P oscillates with an amplitude $h$ and a frequency $\omega$.

We parameterize the pendulum by its length $l$, the amplitude $h$ and
frequency $\omega$ of the driving vertical oscillation of the pivot point, and
the angle $\theta$ of the pendulum with respect to the upright position. Using
$\theta$ as a generalized coordinate, the potential $V_{\mathrm{grav}}=mgy$
where $g=9.81\,\mbox{m}/\mbox{s}^{2}$, and the driving oscillation
parameterized by $\vec{P}(t)=(0,h\sin\omega t)$, the Lagrangian of the system
can be written as

$\displaystyle\mathbf{L}$ $\displaystyle=$ $\displaystyle
ml^{2}\left(\frac{1}{2}\dot{\theta}^{2}-\frac{g}{l}\cos\left(\theta\right)-\frac{h}{l}\dot{\theta}\omega\sin\left(\theta\right)\cos\left(\omega
t\right)\right)$ (1) $\displaystyle\equiv$ $\displaystyle
ml^{2}\,L(\theta,\dot{\theta},t).$

In the second line and throughout the rest of this work we rescale quantities
with dimension of energy (e.g. Lagrangians, potentials, energies) by $ml^{2}$
for notational convenience. The factor $\cos\omega t$ encodes the external
driving oscillation from the pivot point of the pendulum. The range of these
parameters is chosen such that the driving oscillation is a fast vibration of
the pivot with small amplitude, such that $\sqrt{g/l}/\omega\ll 1$ and $h/l\ll
1$. We then rewrite the Lagrangian such that these ratios are described by a
single small parameter $\lambda$, given by the ratio between the natural
frequency of the unperturbed pendulum $\omega_{0}=\sqrt{g/l}$ and the driving
frequency $\omega$. The ratio $h/l$ is written as the same small parameter
$\lambda$ times an $\mathcal{O}(1)$ number $\rho$, i.e.

$\lambda=\frac{1}{\omega}\sqrt{\frac{g}{l}}\,,\qquad\frac{h}{l}=\lambda\rho\,.$
(2)

With these definitions, the Lagrangian can be written as

$\displaystyle
L=\frac{1}{2}\dot{\theta}^{2}-\lambda^{2}\omega^{2}\cos(\theta)-\lambda\omega\rho\dot{\theta}\sin\theta\,\cos\omega
t\,,$ (3)

from which the equation of motion can easily be derived:

$\displaystyle\ddot{\theta}=\lambda\omega^{2}\left(\lambda-\rho\sin\omega
t\right)\sin\theta\,.$ (4)

We now want to study the motion of the pendulum across time scales larger than
the period of the fast oscillation, imagining a situation where the fast
oscillations are too rapid to be resolved by the measurement device and only a
time-averaged motion is seen Kapitza (1951); Stephenson (1908); Landau and
Lifshitz (1960). To this end, one separates the generalized coordinate into a
slowly- and a rapidly-oscillating component $\varphi$ and $\eta$,
respectively,

$\displaystyle\theta(t)=\varphi(t)+\eta(t)\,,$ (5)

where $\varphi(t)$ is assumed to contain frequencies of the order
$\mathcal{O}(\omega_{0})$ and thus only varies slowly over time. It follows
that derivatives acting on the low-frequency mode $\varphi(t)$ are
proportional to the frequency $\omega_{0}=\lambda\omega$ and thus of
$\mathcal{O}(\lambda)$:

$\displaystyle\frac{d^{n}}{dt^{n}}\varphi(t)\sim\mathcal{O}(\omega^{n}\lambda^{n})\varphi(t)\,.$
(6)

On the other hand, the rapid-oscillation mode $\eta(t)$ dominantly contains
frequencies of order $\mathcal{O}(\omega)$. Since the only source of rapid
oscillations is the vibrating pivot, and its amplitude is of
$\mathcal{O}(\lambda)$, the fast oscillating component $\eta(t)$ should itself
be counted as $\mathcal{O}(\lambda)$. Inserting this mode-splitting into the
equation of motion (4), it becomes

$\displaystyle
0=\ddot{\varphi}+\ddot{\eta}-\lambda\omega^{2}\left(\lambda-\rho\sin\omega
t\right)\sin(\varphi+\eta)\,.$ (7)

Using the $\lambda$ power-counting explained above, we expand this expression
in $\lambda$, yielding

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\ddot{\eta}+\lambda\rho\omega^{2}\sin\varphi\sin\omega t$ (8)
$\displaystyle+\,\ddot{\varphi}-\lambda\omega^{2}\left(\lambda\sin\varphi-\rho\,\eta\,\cos\varphi\,\sin\omega
t\right)$ $\displaystyle+\,\mathcal{O}(\lambda^{3})\,,$

where the first line of the right-hand side is of order $\mathcal{O}(\lambda)$
and the second of order $\mathcal{O}(\lambda^{2})$. Because the equation has
to be satisfied at each order in $\lambda$ separately, it splits into two
equations:

$\displaystyle\begin{aligned}
\ddot{\eta}&=-\lambda\rho\omega^{2}\sin\varphi\sin\omega t\,,\\\
\ddot{\varphi}&=\lambda\omega^{2}\left(\lambda\sin\varphi-\rho\,\eta\,\cos\varphi\,\sin\omega
t\right)\,.\end{aligned}$ (9)

In the equation for $\eta$, the slowly varying function $\varphi$ can be
regarded as constant over the time scale of one oscillation of $\eta$. The
equation can thus be integrated directly in $t$, yielding

$\displaystyle\eta=\lambda\rho\sin\varphi\,\sin\omega t\,.$ (10)

The integration constant must be set to zero, since it would constitute a
slowly varying component that is absent from $\eta$ by construction. Inserting
the solution for $\eta$ into the second equation of (9) one finds

$\displaystyle\ddot{\varphi}=\lambda^{2}\omega^{2}\left(1-\rho^{2}\cos\varphi\,\sin^{2}\omega
t\right)\sin\varphi\,.$ (11)

Since the function $\varphi(t)$ varies only slowly with time compared to the
driving oscillation, we rewrite $\sin^{2}\omega t=\frac{1}{2}(1-\cos 2\omega
t)$ and keep only the constant term. This can be interpreted as averaging the
equation over one period of the driving oscillation. The equation of motion
becomes

$\displaystyle\ddot{\varphi}$ $\displaystyle=$
$\displaystyle\lambda^{2}\omega^{2}\left(1-\frac{\rho^{2}}{2}\cos\varphi\right)\sin\varphi$
(12) $\displaystyle\equiv$
$\displaystyle-\frac{d}{d\varphi}V_{\mathrm{eff}}(\varphi)\,,$

which is known as the Kapitza equation. The motion of the pendulum, when
averaged over the rapid oscillations, can be described as the motion in a
time-independent effective potential, which by integrating the right-hand side
of the first line is found to be

$V_{\mathrm{eff}}(\varphi)=\lambda^{2}\omega^{2}\left(\cos\varphi-1+\frac{\rho^{2}}{4}\sin^{2}\varphi\right).$
(13)

The integration constant is chosen such that $V_{\mathrm{eff}}(0)=0$. The
effective potential has a stable equilibrium in the upright position
$\varphi=0$ when $\rho>\sqrt{2}$, which follows from requiring the second
derivative at the origin to be positive,
$V_{\mathrm{eff}}^{\prime\prime}(\varphi)>0$. The potential is shown in Fig. 2
for a stable ($\rho>\sqrt{2}$) and an unstable ($\rho<\sqrt{2}$)
configuration.

Figure 2: Effective potential from eq. (13) for a stable configuration
$\rho>\sqrt{2}$ (gray line), and an unstable configuration $\rho<\sqrt{2}$
(red line).

From the above results, one can derive the averaged pendulum trajectory
$\varphi(t)$, shown for the parameters $\lambda=1/10$ and $\rho=2$ in Fig. 3
as the red line. The figure also shows the numerical solution of the full
system (4), in which the rapidly oscillating components can still be seen. One
notices that the trajectory of $\varphi$ approximates the low-frequency
component of $\theta$ well. But after few periods of the slow oscillation, the
trajectories begin to diverge. To obtain a more precise approximation,
corrections beyond the leading order in $\lambda$ need to be included. We
remark that to compare the two trajectories, one needs to account for the fact
that the initial conditions on $\theta$ in the full system and on $\varphi$ in
the effective system are not identical and need to be matched, as will be
discussed in Sec. II.3 in detail.

Figure 3: Numerical solutions of the exact system (gray line) and the leading-
order effective system (red line) for the parameters $\lambda=1/10$ and
$\rho=2$.

It is interesting to note that the effective Lagrangian is that of a closed
system, since it does not explicitly depend on the time $t$, which is not the
case for the full system. This can be understood from the fact that in the
full system, the time-dependence of the Lagrangian arises through the rapid
oscillations: The vibrating pivot periodically increases and decreases the
energy of the system over the time scale of the fast oscillation. Over much
longer time scales, this effect averages out and the system becomes a closed
one.

The above discussion shows that high-frequency fluctuations that may not be
resolved by observation, can nevertheless change qualitatively the low-
frequency oscillatory behavior, here by creating a stable equilibrium in the
upright position. A similar phenomenon arises in quantum field theory, where
short-distance quantum fluctuations can change the ground-state / vacuum of
the system Coleman and Weinberg (1973). This is currently of great interest
since a competition between the quantum fluctuations of the top quark and
those of the electroweak and Higgs boson determines whether the ground-state
in the effective Higgs potential is stable or metastable Buttazzo  _et al._
(2013).

The derivation of the effective potential (13) presented above is
unsatisfactory for its intransparency as to how it can be systematically
generalized to higher orders in the approximations employed. In both equations
of motion (9) we have made approximations motivated by heuristic scaling
arguments: In the equation for the rapid mode $\eta$, the slow mode $\varphi$
has been treated as static, supported by the scaling argument (6). At higher
orders in $\lambda$, derivatives of $\varphi$ should contribute but their
treatment is unclear. The same issue arises with the averaging procedure
leading to (12).

Previous attempts to solve these issues focused on the improvement of the
time-averaging approach Rahav  _et al._ (2003a); Grundy (2019); Maggia  _et
al._ (2020) often using multi-scale methods Kuehn (2015), and mathematical
analysis of the exact equations of motion Sanz-Serna (2009). Furthermore, many
authors have discussed more accurate approximations of the stability
boundaries Butikov (2011, 2018); Blackburn  _et al._ (1992); Mond  _et al._
(1993); Insperger and Stépán (2002, 2003), but to our knowledge a method to
obtain the systematic expansion of the effective potential in $\lambda$ has
not been discussed before.

In order to address this problem, we draw inspiration from the methodology of
_effective (quantum) field theories_. Throughout all areas of physics, an
important concept in approaching complex systems is to identify the degrees of
freedom that are of subleading importance for the observable of study. For
example, the motion of an electrically charged object in the field generated
by the distribution of charged point particles contained in a finite volume
can be described by the multipole moments of the distribution when the object
is far away compared to the size of the system Jackson (2021). The calculation
of the multipole moments does require knowledge of the “microscopic”
distribution of the point particles, but the details are not relevant for the
object moving on “macroscopic” scales.

The association of effects and their origins with characteristic energy or
length scales is key to the concept of effective field theory. One then
establishes a parametric limit in which the expressions describing the system
can be expanded. The resulting effective theory is then a simpler albeit
approximate description of the system, in which the subleading effects can in
principle be included as power-corrections in an expansion around the
parametric limit to improve the approximation to the desired accuracy.

Effective theories play a prominent role in quantum field theories, where
quantum fluctuations are always present on all scales. To focus on the scales
directly probed by observations of the physical system, one removes
(“integrates out”) the shorter-ranged quantum fluctuations to arrive at an
effective description in terms of the relevant degrees of freedom Weinberg
(1979); Polchinski (1992); Manohar (2020). The inverted pendulum with
$\omega\gg\omega_{0}$ provides a classical analogy of removing a small time
scale, provided by the inverse of the rapid oscillation frequency of the pivot
point. The condition $h\ll l$ implies that the rapid oscillations are “weak”,
so that they can be treated perturbatively in the small parameter
$\lambda=\omega_{0}/\omega$.

An important aspect in finding effective descriptions is the ability to go
beyond the leading approximations in a systematic fashion by including higher-
order terms in the power-expansion in $\lambda$. While simple from a
conceptual point of view, the treatment of subleading-power effects often
brings technical subtleties not present at the leading power of the expansion.

In the following, we provide a derivation of the effective Lagrangian for the
inverted pendulum that differs from the approach demonstrated above, is close
in spirit to the mode decomposition common to effective quantum field
theories, and can be systematically extended to any order in the expansion
parameter $\lambda$. We first show the derivation of the first non-trivial
power correction in detail and then give results based on a higher-power
analysis.

##  II Subleading-Power Analysis

###  II.1 Mode Expansion

We now present a more systematic approach to the problem. Our ansatz is to
make the frequency-separation manifest by factoring out all fast oscillations,
which come in multiples of the large frequency $\omega$, from the generalized
coordinate in terms of the modulated Fourier expansion Sanz-Serna (2009)

$\displaystyle\theta(t)=\sum_{n=0}^{\infty}\left[\Psi_{n}(t)\sin(n\omega
t)+\Phi_{n}(t)\cos(n\omega t)\right]\,.$ (14)

By construction, the coefficient mode functions $\Psi_{n}$ and $\Phi_{n}$ of
the rapidly varying trigonometric functions only vary slowly with the time
$t$, specifically with time scales of the order $\mathcal{O}(1/\omega_{0})$.
The low-frequency mode $\varphi$ from the previous section is then the
contribution from $n=0$. Including only discrete values of $n$ is justified by
the fact that the sole origin of high-frequency oscillations is the sinusoidal
source term proportional to $\lambda\cos\omega t$ in the Lagrangian (3). We
now note that the source term coupling the coordinate $\theta$ to the high-
frequency external vibration comes with one power of $\lambda$, implying that
the amplitude of $n$-th harmonic of $\omega$ is naturally suppressed by
$\lambda^{n}$. We also deduce that the even-$n$ terms can only involve
$\Phi_{n}$ whereas the odd-$n$ terms only have $\Psi_{n}$ coefficients. To
solve the equations for $\Phi_{n},\Psi_{n}$ successively in powers of
$\lambda$, we expand the mode functions for $n>0$ in the power-counting
parameter,

$\displaystyle\Psi_{n}(t)$ $\displaystyle=\sum_{m=0}^{\infty}\psi_{nm}(t)\,,$
$\displaystyle\Phi_{n}(t)$
$\displaystyle=\sum_{m=0}^{\infty}\varphi_{nm}(t)\,,$ (15)

and identify the zero mode (recall $\Psi_{0}=0$)as

$\Theta(t)\equiv\Phi_{0}(t)\,.$ (16)

In these series, the $m$-th terms scale as

$\psi_{nm}(t),\,\varphi_{nm}(t)\ \sim\ \lambda^{n+m}\,,$ (17)

and since by assumption these functions are slowly varying, their $k$-th
derivatives count as

$\displaystyle\psi_{nm}^{(k)}(t),\,\varphi_{nm}^{(k)}(t)\ \sim\
\lambda^{n+m+k}\,.$ (18)

For later convenience, we introduce the abbreviations

$\displaystyle\sin(n\cdot\Theta(t))=\mathrm{s}_{n\Theta}\,,\;\;\cos(n\cdot\Theta(t))=\mathrm{c}_{n\Theta}\,.\quad$
(19)

We emphasize that unlike the fast components of $\theta(t)$, the zero mode
$\Theta(t)$ and consequently $\mathrm{s}_{n\Theta}$, $\mathrm{c}_{n\Theta}$
are not expanded in $\lambda$ in the following.

We now insert eq. (14) into the Lagrangian (3) and expand to the desired order
in $\lambda$. Since all factors of $\sin\omega t$ and $\cos\omega t$ are
explicit in the mode expansion of the coordinate $\theta$, the Lagrangian
itself automatically decomposes into Fourier modes, casting it into the form:

$\displaystyle L$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\left[A_{k}(\Phi,\dot{\Phi},\Psi,\dot{\Psi})\cdot\cos(k\omega
t)\right.$ (20) $\displaystyle\left.\hskip
17.07182pt+\,B_{k}(\Phi,\dot{\Phi},\Psi,\dot{\Psi})\cdot\sin(k\omega
t)\right].$

The coefficients of the rapid oscillations $\cos(k\omega t)$ and $\sin(k\omega
t)$ in this equation are not to be confused with the ones in the mode
decomposition of the generalized coordinate $\theta$ appearing in the
decomposition (14), as the full Lagrangian (3) is non-linear in $\theta$.
Products of different fast Fourier modes of $\theta$ then decompose into a sum
of terms in the Lagrangian (20) after making use of standard trigonometric
identities.

The Lagrangian (20) is still equivalent to the one of the full system, and
describes its dynamics at all frequency scales. The slowly oscillating
behavior of the system, which we are ultimately interested in, is given by the
$k=0$ term in the sum, $L_{k=0}\equiv A_{0}$. This expression contains
contributions from the $n=0$ modes of $\theta$ as well as terms in which the
higher modes combine to yield a slowly-varying function. This is equivalent to
the averaging procedure in the previous section, where we expanded the
solution of the low-frequency mode (11) using trigonometric identities and
kept only the term without $\cos(k\omega t)$. To determine the mode functions
$\psi_{nm}(t)$, $\varphi_{nm}(t)$, we insert the mode expansion into the
equation of motion of the full system, and solve it separately mode by mode
and order by order in the expansion in $\lambda$.

We start at the equation of motion (4) of the full system,

$\displaystyle\ddot{\theta}=\lambda^{2}\omega^{2}\sin\theta-\lambda\rho\omega^{2}\sin\theta\sin\omega
t\,,$ (21)

and insert eqs. (14) and (15) into the above expression. We then expand in
powers of $\lambda$ according to the defined scaling rules (17), (18) for all
the objects. For example, up to order $\lambda^{2}$, the equation of motion
becomes:

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\lambda\omega^{2}\sin(\omega
t)\left[\rho\mathrm{s}_{{}\Theta}-\frac{\psi_{10}}{\lambda}\right]$
$\displaystyle+\,\lambda^{2}\left(\left[\frac{\ddot{\Theta}}{\lambda^{2}}-\omega^{2}\mathrm{s}_{{}\Theta}+\frac{\rho^{2}\omega^{2}}{2}\mathrm{c}_{{}\Theta}\frac{\psi_{10}}{\lambda}\right]\right.$
$\displaystyle+\left.\omega\cos(\omega
t)\left[\frac{2\dot{\psi}_{10}}{\lambda^{2}}-\omega\frac{\varphi_{11}}{\lambda^{2}}\right]\right.$
$\displaystyle-\left.\frac{\omega^{2}}{2}\cos(2\omega
t)\left[\rho\mathrm{c}_{{}\Theta}\frac{\psi_{10}}{\lambda}+\frac{8\varphi_{20}}{\lambda^{2}}\right]\right)+\mathcal{O}(\lambda^{3})\,.$

The equation needs to be fulfilled at each order in $\lambda$ separately, as
well as for each harmonic of $\omega$. In eq. (II.1) this means that each
square bracket has to vanish independently, defining a system of equations.
One then solves them at each order in $\lambda$, inserting the solutions in
the equations at the next higher order until all mode functions $\psi_{nm}$,
$\varphi_{nm}$ are determined. In the above example, this means we start by
setting the expression in square brackets in the first line to zero, yielding

$\displaystyle\psi_{10}=\lambda\rho\mathrm{s}_{{}\Theta}\,.$ (23)

Inserting this solution into the $\mathcal{O}(\lambda^{2})$ terms of eq.
(II.1) and requiring each bracket to vanish, one finds the solution for the
next-to-leading order modes:

$\displaystyle\begin{aligned}
\varphi_{20}(t)&=-\frac{\lambda^{2}\rho^{2}}{8}\mathrm{c}_{{}\Theta}\mathrm{s}_{{}\Theta}\,,\\\
\varphi_{11}(t)&=\frac{2\lambda\rho}{\omega}\dot{\Theta}\mathrm{c}_{{}\Theta}\,.\end{aligned}$
(24)

The square bracket in the second line of eq. (II.1) contains no further
information about any $\varphi_{nm}$, $\psi_{nm}$ and can be discarded – it
contains instead the leading-power equation of motion of the zero-mode
$\Theta(t)$. Using the above results, we can insert them into $L_{k=0}$ to
obtain the _effective Lagrangian_ up to order $\mathcal{O}(\lambda^{4})$:

$\displaystyle\begin{aligned}
L_{\mathrm{eff}}=&\frac{1}{2}\dot{\Theta}^{2}\left(1-\frac{3}{2}\lambda^{2}\rho^{2}\mathrm{c}_{{}\Theta}^{2}\right)-\lambda^{2}\omega^{2}\left(\mathrm{c}_{{}\Theta}-\frac{\rho^{2}}{8}\mathrm{c}_{2\Theta}\right)+\frac{\lambda^{4}\omega^{2}\rho^{2}}{16}\left(\mathrm{c}_{{}\Theta}-\frac{\rho^{2}}{2}\mathrm{c}_{2\Theta}-\mathrm{c}_{3\Theta}+\frac{5}{32}\rho^{2}\mathrm{c}_{4\Theta}\right)\\\
&+\lambda^{2}\rho^{2}\mathrm{s}_{{}\Theta}\ddot{\Theta}\left(\frac{5}{64}\lambda^{2}\rho^{2}\mathrm{c}_{3\Theta}-\lambda^{2}\mathrm{c}_{{}\Theta}^{2}-\frac{5}{\omega^{2}}\mathrm{c}_{{}\Theta}\dot{\Theta}^{2}\right)+\frac{\lambda^{2}\rho^{2}}{4\omega^{2}}\mathrm{c}_{{}\Theta}^{2}\ddot{\Theta}^{2}+\frac{3\lambda^{2}\rho^{2}}{2\omega^{2}}\mathrm{c}_{{}\Theta}^{2}\dot{\Theta}\Theta^{(3)}+\ldots\end{aligned}$
(25)

The second line adds a subset of $\mathcal{O}(\lambda^{6})$ terms to be
discussed below. This effective Lagrangian determines the dynamics of the
system when the fast oscillations at the scale $\omega$ are not resolved by
the measurement. By determining the equation of motion for $\Theta$ and
solving it, we obtain the motion of the pendulum when the fast fluctuations
are averaged over. Some comments are in order.

First, the above procedure can be generalized to an arbitrary order in
$\lambda$, and at any order one finds the equations determining the mode
functions to be purely algebraic and never differential. This is due to the
fact that derivatives of mode functions are always power-suppressed. Thus,
equations involving derivatives of a given mode function appear at higher
orders in $\lambda$ than equations involving the same function without
derivatives, but the corresponding mode is always already determined at lower
orders through an algebraic equation.

Second, the effective Lagrangian does not explicitly depend on time, and the
averaged system remains conservative to any order in the $\lambda$ expansion.
This is not true for the full system, as the pivot oscillation periodically
injects energy to and removes energy from the system. However, within the
domain of validity of the expansion in $\lambda$, the effect averages out and
rapid oscillations do not transfer energy to the pendulum. The system is
effectively closed.

Third, at $\mathcal{O}(\lambda^{4})$ the effective Lagrangian contains a
coordinate-dependent kinetic term of the form
$\dot{\Theta}^{2}\mathrm{c}_{{}\Theta}^{2}$, which is a general feature of the
slow-motion dynamics in rapidly oscillating potentials Rahav  _et al._ (2003a,
b). Furthermore, starting at order $\mathcal{O}(\lambda^{6})$, the effective
Lagrangian contains higher derivatives of $\Theta$. These terms are displayed
explicitly in the second line of (25). In general, at
$\mathcal{O}(\lambda^{2n+2})$ ($n>0$) any monomial of up to the $(2n-1)$-th
derivative of $\Theta$ with a total of $2n$ derivatives can appear. For
example, at order $\mathcal{O}(\lambda^{6})$:

$\dot{\Theta}^{4},\dot{\Theta}^{2}\ddot{\Theta},\dot{\Theta}\Theta^{(3)},\ddot{\Theta}^{2}.$
(26)

Ordinarily, the presence of such terms in the Lagrangian indicates an
Ostrogradsky instability of the system, i.e. that the system’s energy is
unbounded from below Ostrogradsky (1850). However, in an effective description
higher-derivative terms do not lead to an inconsistency if they are power-
suppressed. As such they must be treated as perturbations to the leading-power
Lagrangian that does not involve higher derivatives. Using power-counting
arguments, we will show below that for the inverted pendulum, the higher-
derivative terms and higher powers of $\dot{\Theta}^{2}$ can be traded for
subleading terms with fewer derivatives by means of the equations of motion
and energy conservation.

The emergence of higher-derivative terms in the systematic treatment of the
slow-motion dynamics of the inverted pendulum presents another similarity to
the case of effective relativistic quantum field theories, in which the
effective Lagrangian can contain kinetic terms with more derivatives than the
canonical one, $\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi$. For example,
a four-derivative term
$(\partial_{\mu}\phi\partial^{\mu}\phi)^{2}/\Lambda^{2}$, where $\Lambda$
represents the large momentum scale, would cause an instability and further
imply additional tachyonic modes violating relativistic causality. However,
for slowly varying fields compared to the scale $\Lambda$, the power-
suppression of this term implies that it should be eliminated by perturbative
field redefinitions. At leading order in the power-counting of the theory,
this procedure is equivalent to using the equations of motion of the field.

These remarks also clarify a potentially troublesome point, namely that the
presence of higher-derivatives of $\Theta$ in the effective Lagrangian for the
inverted pendulum does not imply that the solution to the effective equation
of motion needs initial conditions beyond those for $\Theta$ and
$\dot{\Theta}$. Instead, the initial conditions follow from matching to those
of the full system, as detailed below.

###  II.2 Effective Potential

We now show that the higher derivatives and higher powers of
$\dot{\Theta}^{2}$ can be eliminated, and the effective Lagrangian converted
into the canonical form

$\displaystyle
L_{\mathrm{eff}}=\frac{1}{2}\dot{\Theta}^{2}-V_{\mathrm{eff}}(\Theta,E_{\rm
tot})\,$ (27)

with an effective potential, which beyond the well-known leading approximation
(13) turns out to depend on the total energy of the system.

Since the effective Lagrangian in the previous subsection does not depend
explicitly on time, the Noether charge associated with the time-translation
symmetry defines the conserved energy. For a Lagrangian depending on higher
derivatives of the generalized coordinate $\Theta$, the relevant expression is
Anderson (1973)

$\displaystyle E_{\rm
tot}\\!=\sum_{k=1}^{\infty}\sum_{j=0}^{k-1}(-1)^{j}\Theta^{(k-j)}\frac{d^{j}}{dt^{j}}\frac{\partial
L_{\mathrm{eff}}}{\partial\Theta^{(k)}}\\!-L_{\mathrm{eff}}\,,$ (28)

which generalizes the standard expression for the total energy. From the
effective Lagrangian (25), we obtain

$\displaystyle\begin{aligned}
E_{\mathrm{tot}}=&\frac{1}{2}\dot{\Theta}^{2}\left(1-\frac{3\lambda^{2}\rho^{2}}{2}\mathrm{c}_{{}\Theta}^{2}\right)+\lambda^{2}\omega^{2}\left(\mathrm{c}_{{}\Theta}-\frac{\rho^{2}}{8}\mathrm{c}_{2\Theta}\right)+\frac{\lambda^{4}\rho^{2}\omega^{2}}{16}\left(-\mathrm{c}_{{}\Theta}+\frac{\rho^{2}}{2}\mathrm{c}_{2\Theta}+\mathrm{c}_{3\Theta}-\frac{5\rho^{2}}{32}\mathrm{c}_{4\Theta}\right)\\\
&-\frac{5\lambda^{2}\rho^{2}}{4\omega^{2}}\mathrm{c}_{{}\Theta}^{2}\ddot{\Theta}^{2}+\frac{5\lambda^{2}\rho^{2}}{2\omega^{2}}\mathrm{c}_{{}\Theta}^{2}\dot{\Theta}\Theta^{(3)}-\frac{5\lambda^{2}\rho^{2}}{\omega^{2}}\dot{\Theta}^{2}\mathrm{c}_{{}\Theta}\mathrm{s}_{{}\Theta}\ddot{\Theta}+\ldots\,,\end{aligned}$
(29)

where in the second line again only the higher derivative terms at
$\mathcal{O}(\lambda^{6})$ are shown. We note that the leading-order
expression for the energy does not depend on higher derivatives and higher
powers of $\dot{\Theta}^{2}$. It is therefore possible to define an effective
potential by successively solving

$\displaystyle E_{\rm
tot}=\frac{1}{2}\dot{\Theta}^{2}+V_{\mathrm{eff}}(\Theta,E_{\rm tot})$ (30)

order by order in $\lambda$ to eliminate the higher derivatives of $\Theta$
and powers of $\dot{\Theta}^{2}$, which yields111 The superscript $(n)$ on
$V_{\rm eff}$, $E_{\rm tot}$ and $L_{\rm eff}$ here and below denotes the
$\mathcal{O}(\lambda^{n})$ term in the corresponding quantity.

$\displaystyle V_{\mathrm{eff}}(\Theta,E_{\rm
tot})\equiv\sum_{n=1}^{\infty}V_{\mathrm{eff}}^{(2n)}(\Theta,E_{\rm tot})\,.$
(31)

The effective potential obtained through this procedure depends on the total
energy of the effective system. Expanding $E_{\rm tot}$ and $V_{\mathrm{eff}}$
in powers of $\lambda$, $V_{\mathrm{eff}}^{(2n)}$ will generally depend on
$E_{\mathrm{tot}}^{(2m)}$ with $m<n$.

Here we demonstrate the procedure up to the next-to-leading order in the
$\lambda$ expansion. Beginning with the leading order, eq. (28) becomes

$\displaystyle\begin{aligned} E_{\rm tot}^{(2)}&=\dot{\Theta}\frac{\partial
L^{(2)}_{\mathrm{eff}}}{\partial\dot{\Theta}}-L^{(2)}_{\mathrm{eff}}\,.\end{aligned}$
(32)

Combining this with eq. (30), one obtains

$\displaystyle
V^{(2)}_{\mathrm{eff}}=\lambda^{2}\omega^{2}\left(\mathrm{c}_{{}\Theta}+\frac{\rho^{2}}{4}\mathrm{s}_{{}\Theta}^{2}\right)+\delta
V^{(2)}_{\mathrm{eff}}\,,$ (33)

where $\delta V_{\mathrm{eff}}^{(2)}=-\lambda^{2}\omega^{2}$ is a constant
shift introduced in order to adjust the potential such that
$V_{\mathrm{eff}}(0)=0$. This order reproduces the familiar result (13).

At the next order, from the $\mathcal{O}(\lambda^{4})$ terms of (29), the
procedure yields the following equation to solve

$\displaystyle V_{\mathrm{eff}}^{(4)}+\delta
V_{\mathrm{eff}}^{(4)}\stackrel{{\scriptstyle!}}{{=}}-\frac{3\lambda^{2}\rho^{2}}{4}\mathrm{c}_{{}\Theta}^{2}\dot{\Theta}^{2}$
$\displaystyle+\,\frac{\lambda^{4}\rho^{2}\omega^{2}}{16}\,\bigg{(}\\!\\!-\mathrm{c}_{{}\Theta}+\frac{\rho^{2}}{2}\mathrm{c}_{2\Theta}+\mathrm{c}_{3\Theta}-\frac{5\rho^{2}}{32}\mathrm{c}_{4\Theta}\\!\bigg{)}.\;\quad$
(34)

The $\dot{\Theta}^{2}$ term is eliminated using eq. (30) at
$\mathcal{O}(\lambda^{2})$, expressing it as $\dot{\Theta}^{2}=2\,(E_{\rm
tot}^{(2)}-V_{\mathrm{eff}}^{(2)})$. Tweaking the offset $\delta
V_{\mathrm{eff}}^{(4)}$ such that $V_{\rm eff}(0)=0$, we obtain

$\displaystyle V^{(4)}_{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle\lambda^{4}\rho^{2}\omega^{2}\left\\{\frac{17}{16}\mathrm{c}_{{}\Theta}-\frac{3}{4}\mathrm{c}_{2\Theta}+\frac{7}{16}\mathrm{c}_{3\Theta}\right.$
(35)
$\displaystyle\left.-\,\frac{3}{4}+\rho^{2}\left(\frac{1}{32}\mathrm{c}_{2\Theta}-\frac{29}{512}\mathrm{c}_{4\Theta}+\frac{13}{512}\right)\right\\}$
$\displaystyle+\,\frac{3}{2}\lambda^{2}\rho^{2}\mathrm{s}_{{}\Theta}^{2}E_{\mathrm{tot}}^{(2)}\,.$

When working beyond the fourth order in $\lambda$, one has to eliminate
$\dot{\Theta}^{2}$ consistently to sufficiently high order. Specifically, to
obtain the $\mathcal{O}(\lambda^{2n})$ term, $V_{\mathrm{eff}}^{(2n)}$, in the
effective potential, one subtracts $\frac{1}{2}\dot{\Theta}^{2}$ from the
right-hand side of (29) and eliminates powers of $\dot{\Theta}^{2}$ in the
remaining expression by the solution of (30) to order $\lambda^{2n-2}$, which
includes the previously determined potential terms $V_{\mathrm{eff}}^{(2k)}$
with $k<n$.

Starting from $\mathcal{O}(\lambda^{6})$, higher derivative terms appear,
which can be likewise eliminated by taking the appropriate number of
derivatives of (30) and proceeding as described above for the higher powers of
$\dot{\Theta}^{2}$.

The total energy, which first appears at $\mathcal{O}(\lambda^{4})$ in the
effective potential is determined by the initial conditions at some time
$t_{i}$, thus expressing it through $\Theta(t_{i})\equiv\Theta_{i}$ and
$\dot{\Theta}(t_{i})\equiv\dot{\Theta}_{i}$. At $\mathcal{O}(\lambda^{4})$,
only

$\displaystyle E_{\rm
tot}^{(2)}=\frac{1}{2}\dot{\Theta}^{2}_{i}+V_{\mathrm{eff}}^{(2)}(\Theta_{i})$
(36)

is needed. The initial condition of the effective system is related to that of
the full one as will be explained in the following section.

With the above results in place, we have an expression for $V_{\mathrm{eff}}$
that defines a genuine potential, i.e. a function of only $\Theta$ and no
derivatives. The result can easily be generalized to an arbitrary order in the
expansion in $\lambda$ by iteratively applying the procedure above, order by
order. We give these results in Appendix A to $\mathcal{O}(\lambda^{8})$.

###  II.3 Initial-Condition Matching

When deriving the pendulum motion $\Theta(t)$ from the effective potential and
comparing it to the solution of the full system, one has to account for the
fact that the initial conditions for the full and the effective system differ,
since the full solution still contains the rapid driving oscillations. For
example, the full system could have no initial velocity, because at the
initial time the slow and fast oscillation accidentally cancel out. The
effective system, containing only the slow mode, would then have a non-zero
initial velocity. Therefore, when attempting to compare solutions, as shown in
Fig. 3, the initial conditions on the full trajectory $\theta(t)$ and the
effective one $\Theta(t)$ need to be carefully related to each other.

As an aside, in applications of effective field theory to scattering in
quantum field theory, the matching of initial conditions is unnecessary, as
the asymptotic in- and out-states of free particles at $t\to\mp\infty$ are
expressed in terms of the same degrees of freedom. The situation discussed
here for the inverted pendulum is analogous to the matching of correlation
functions of fields, when the full-theory and effective theory field variables
are different. This arises, for example, in the effective theory for the long-
wavelength fluctuations of a massless scalar field on the De Sitter space-time
background Cohen and Green (2020); Cohen  _et al._ (2021).

We begin by rewriting the mode expansion (5) and solving it for the effective
coordinate $\Theta$ at a given initial time $t_{i}$,

$\displaystyle\begin{aligned} \Theta_{i}&=\theta_{i}\\\
&-\sum_{n=1}^{\infty}\left[\Psi_{n}(t_{i})\sin(n\omega
t_{i})+\Phi_{n}(t_{i})\cos(n\omega t_{i})\right]\,.\end{aligned}$ (37)

Here we separated the zero mode $\Theta$ from the higher ones and the
coordinates with subscript $i$ are meant to be evaluated at the initial time,
$\Theta_{i}=\Theta(t_{i})$. Inserting the obtained solutions for the various
mode functions $\Psi_{n}$ and $\Phi_{n}$ allows one to relate the initial
conditions $\theta_{i}$ to those of the effective coordinate $\Theta_{i}$.
Taking derivatives does the same for the initial angular velocities:

$\displaystyle\begin{aligned} \dot{\Theta}_{i}&=\dot{\theta}_{i}\\\
&-\sum_{n=1}^{\infty}\Big{[}\dot{\Psi}_{n}(t_{i})\sin(n\omega
t_{i})+\dot{\Phi}_{n}(t_{i})\cos(n\omega t_{i})\\\
&+n\omega\big{(}\Psi_{n}(t_{i})\cos(n\omega t_{i})-\Phi_{n}(t_{i})\sin(n\omega
t_{i})\big{)}\Big{]}\,.\end{aligned}$ (38)

At $\mathcal{O}(\lambda)$, the above results translate to

$\displaystyle\begin{aligned}
\Theta_{i}&=\theta_{i}-\lambda\rho\sin\Theta_{i}\sin\omega t_{i}\,,\\\
\dot{\Theta}_{i}&=\dot{\theta}_{i}-\lambda\rho\omega\sin\Theta_{i}\cos\omega
t_{i}\,.\end{aligned}$ (39)

Note that the second equation does not come with a term involving
$\dot{\Theta}_{i}$, as $\dot{\Theta}_{i}$ is suppressed by an additional power
of $\lambda$. Both sides of the equations contain the initial conditions
$\Theta_{i}$, $\dot{\Theta}_{i}$ which we aim to solve for. However, they
always appear on the right-hand sides in power-suppressed terms. Therefore,
the system can be solved iteratively by expanding both $\Theta_{i}$ and
$\dot{\Theta}_{i}$ in powers of $\lambda$. At $\mathcal{O}(\lambda)$, this
amounts to setting

$\displaystyle\begin{aligned}
\Theta_{i}&=\theta_{i}-\lambda\rho\sin\theta_{i}\sin\omega t_{i}\,,\\\
\dot{\Theta}_{i}&=\dot{\theta}_{i}-\lambda\rho\left(\omega\sin\theta_{i}\cos\omega
t_{i}\right)\,.\end{aligned}$ (40)

We give the results to higher order in Appendix B. It is important to recall
that, at higher orders in power-counting, higher derivatives are introduced
again. They can be eliminated with the method described in the last section.
With these results, one can derive numerical solutions for the effective
system and compare them directly to the ones for the full system, as shown in
Fig. 3 above and Fig. 5 below.

The initial conditions are also needed to eliminate the $E_{\rm tot}$ terms in
the effective potential, as mentioned in the previous section for the case of
$E_{\rm tot}^{(2)}$ in $V_{\mathrm{eff}}^{(4)}$, where one only needed the
leading-power expression (36). In general, at $\mathcal{O}(\lambda^{2n})$, eq.
(30) evaluated at initial time $t_{i}$ gives

$\displaystyle\dot{\Theta}_{i}^{2}=2\sum_{j=1}^{n}\left(E_{\rm
tot}^{(2j)}-V_{\mathrm{eff}}^{(2j)}(\Theta_{i},E_{\rm tot})\right)\,,$ (41)

where the argument of $V_{\mathrm{eff}}^{(2j)}$ involves the expansion terms
$E_{\mathrm{tot}}^{(2m)}$ with $m$ taking values up to $j-1$. Assuming that
$E_{\mathrm{tot}}^{(2m)}$ is known up to $m=n-1$, the above equation can be
solved for $E_{\rm tot}^{(2n)}(\Theta_{i},\dot{\Theta}_{i})$. Finally,
$\Theta_{i},\dot{\Theta}_{i}$ are expressed in terms of the initial conditions
$\theta_{i},\dot{\theta}_{i}$ of the full system to order $\lambda^{2n}$ as
described above to find $E_{\rm tot}^{(2n)}(\theta_{i},\dot{\theta}_{i})$
after re-expanding the entire expression in $\lambda$.

##  III Discussion of Results

Figure 4: Effective potentials $V_{\mathrm{eff}}^{NnLO}$, where $n$ denotes
the order in power-counting, with $\lambda=1/4$ and $\rho=2$ and initial
conditions $\theta(0)=0.05$, $\dot{\theta}(0)=-0.03\mathrm{s}^{-1}$.

###  III.1 Effective potential and low-frequency motion

We begin by comparing the effective potentials at different orders up to order
$\mathcal{O}(\lambda^{8})$. The explicit expressions are provided in App. A.
To display the result, the dependence of the effective potential on the total
energy $E_{\rm tot}^{(n)}$ is eliminated through initial-condition matching.
For the purpose of illustration we choose the initial conditions
$\theta(0)=0.05$, $\dot{\theta}(0)=-0.03\mathrm{s}^{-1}$, and $\lambda=1/4$,
$\rho=2$. Fig. 4 shows the effective potentials $V_{\rm eff}^{\rm LO}$,
$V_{\rm eff}^{\rm NLO}$, $V_{\rm eff}^{\rm N3LO}$ obtained including the
correction to $\mathcal{O}(\lambda^{2})$, $\mathcal{O}(\lambda^{4})$ and
$\mathcal{O}(\lambda^{8})$, respectively, drawn in gray, blue and red,
respectively. We see that the leading-power result deviates noticeably from
the higher-order ones, which are already close to each other, indicating good
convergence of the expansion for the chosen parameters.

Figure 5: Motion of the pendulum from the subleading-power analysis, using the
same parameters as for Fig. 4.

From the effective potential we can next derive the effective motion of the
pendulum and compare it against the one obtained from the solution of eq. (4),
which still contains the rapid oscillations. To do so, the initial conditions
must be carefully adjusted order by order in $\lambda$ as detailed in Sec.
II.3. We show the averaged motions in Fig. 5 at leading order, $\Theta^{\rm
LO}(t)$, in red and $\Theta^{\rm N3LO}(t)$, obtained from the
$\mathcal{O}(\lambda^{8})$-accurate potential $V_{\rm eff}^{\rm N3LO}$, in
blue. For comparison, we show the full solution $\theta(t)$ and its time-
average $\bar{\theta}_{2\omega_{0}}(t)$, defined as $\theta(t)$ after applying
a low-pass filter with cut-off frequency $2\omega_{0}$ as the gray dashed and
solid lines, respectively. In order to quantify the deviation of the effective
solution from the exact one, we also determine the averaged square difference,

$\displaystyle\Delta^{(\mathrm{N}n\mathrm{LO})}=\frac{1}{T}\int_{0}^{T}dt\,(\bar{\theta}_{2\omega_{0}}(t)-\Theta^{(\mathrm{N}n\mathrm{LO})}(t))^{2}\,,$
(42)

and show the result in Fig. 6, also including different values for the
expansion parameter $\lambda$, but keeping $\rho$ and the initial condition
values unchanged. Rapid improvement with increasing order of approximation is
seen in all three cases.

Figure 6: Deviations $\Delta^{(\mathrm{N}n\mathrm{LO})}$ of the effective and
full solutions shown in Fig. 5, as defined in eq. (42). For comparison, points
for other values of $\lambda$ are also included.

###  III.2 Stability of the upright equilibrium

Stability of the upright equilibrium requires

$\displaystyle\left.\frac{d^{2}}{d\Theta^{2}}V_{\mathrm{eff}}(\Theta,E_{\rm
tot})\right|_{\Theta=0}>0\,.$ (43)

From the above expressions for $V_{\mathrm{eff}}^{(2)}$,
$V_{\mathrm{eff}}^{(4)}$, we obtain

$\displaystyle\rho^{2}>2+\frac{7}{4}\lambda^{2}-\frac{12}{\omega^{2}}E_{\rm
tot}^{(2)}+\mathcal{O}(\lambda^{4})\,,$ (44)

which includes the $\mathcal{O}(\lambda^{2})$ correction to $\rho^{2}>2$.
Interestingly, the inverted pendulum becomes stable already for smaller
amplitude $h$ of the pivot point oscillation when the energy of the system
increases. Yet, since the energy dependence is a subleading-power effect in
$\lambda$, the depth of the potential well does not increase significantly.

If we are content with studying the frequency of a very small oscillation
around the upright equilibrium point, we can obtain results to higher orders
in $\lambda$.222The equation of motion of the full system in this limit takes
the form of the well-known Mathieu equation. For simplicity, we set $E_{\rm
tot}=0$ in the effective potential. In the limit of small amplitudes, the
effective potential then directly yields the frequency of the (now sinusoidal)
oscillation around the upright equilibrium, using

$\displaystyle
V_{\mathrm{eff}}=\frac{\omega_{\mathrm{eff}}^{2}\Theta^{2}}{2}+\mathcal{O}(\Theta^{3})\,.$
(45)

From this we obtain the effective frequency $\omega_{\mathrm{eff}}$ of small-
angle oscillations, which is expanded into a power-series in $\lambda$ as

$\displaystyle\omega_{\mathrm{eff}}^{2}=\omega^{2}\sum_{n=1}^{\infty}\lambda^{2n}c_{\omega}^{(2n)}\,.$
(46)

Using computer algebra software, we calculated the expansion of
$\omega_{\mathrm{eff}}^{2}$ up to order $\mathcal{O}(\lambda^{28})$, which
corresponds to thirteen additional terms of the series beyond the leading
approximation. The first three coefficients are found to be

$\displaystyle\begin{aligned} c_{\omega}^{(2)}&=\frac{\rho^{2}}{2}-1\,,\\\
c_{\omega}^{(4)}&=\frac{25}{32}\rho^{4}-2\rho^{2}\,,\\\
c_{\omega}^{(6)}&=\frac{1169}{576}\rho^{6}-\frac{273}{32}\rho^{4}+8\rho^{2}\,.\end{aligned}$
(47)

The complete set of calculated coefficients is given in App. C. Their
magnitude relative to the leading term is shown in Fig. 7. The individual
terms are seen to fall off roughly exponentially with growing $n$.

Figure 7: Relative corrections to the effective frequency as defined in eq.
(46) with $\rho=2$.

##  IV Convergence Behavior

Equipped with a method to compute an effective Lagrangian to (in principle)
any order in $\lambda$, it is interesting to investigate whether or not the
series expansions converge. We consider the small oscillation frequency
$\omega_{\rm eff}$, for which the first 14 terms of the series are at our
disposal. In quantum mechanics (see Bender and Wu (1969) for the anharmonic
oscillator) and quantum field theory Dyson (1952), perturbation expansions for
energy eigenvalues or correlation functions usually have zero radius of
convergence and the series are at best asymptotic. Inspection of the series
coefficients for the classical mechanics problem at hand, however, suggests
that there is a finite radius of convergence in $\lambda$. Failure of an
expansion to converge often indicates that a qualitatively new physical
phenomenon emerges at the convergence boundary. Indeed, one must expect the
$\lambda$ expansion to fail for too large $\lambda=\omega_{0}/\omega$ or
$\lambda\rho=h/l$, and the effective description to break down, as the
inverted pendulum exhibits parametric resonance Butikov (2011, 2018); Mond
_et al._ (1993); Blackburn  _et al._ (1992); Insperger and Stépán (2002,
2003), in which case the motion of the pendulum can certainly not be described
by a time-independent effective Lagrangian with a conserved energy. This
happens when the natural frequency $\omega_{0}$ of the pendulum and the
driving frequency $\omega$ are close to each other, in which case there is no
frequency separation and energy from the driver is transferred to the pendulum
until it escapes the potential well around the upright equilibrium.

To establish a baseline to which to compare the results from the effective
description to high order in $\lambda$, we first investigate the stability of
the full system by conducting a parameter scan of the trajectories in the
$(\lambda\rho,\rho)$-plane. To simplify the following discussion and to make
it more comparable to existing analyses on the Mathieu equation, we will focus
solely on the small angle limit from now on. At each point, the full system is
solved numerically with initial values $\theta(0)=0$ and
$\dot{\theta}(0)=10^{-5}/\mathrm{s}$, and the time until the oscillation has
amplified such that the pendulum tips over from the upright position (meaning
it passes the point $|\theta|>\pi/2$) is measured. We solve the system until
$t_{\mathrm{max}}=300\mathrm{s}$ and take $\omega$ for each parameter point
such that $\omega_{0}=\sqrt{g/l}=3.14/\mathrm{s}$ (recall
$\omega=\omega_{0}/\lambda$). If it does not tip over within the time scale of
simulation, the point is considered stable. Choosing $t_{\rm
max}=300\mathrm{s}$ turned out to be a reasonable criterion, because the
typical tip-over time was orders of magnitude smaller than $t_{\rm max}$,
apart from a region narrower than the provided resolution.

We show the results of this simulation in Fig. 8, where the coloring indicates
the time at which the system tips over. Brighter colors denote solutions that
tip over more quickly. The black wedge-shaped region is the parameter space in
which the system is considered stable by the means of our simulation. The tip-
over time exhibits interesting patterns in the instability region, which,
however, are not our concern here. Instead, in what follows, we show how the
effective description can be used to determine the boundary between the stable
and unstable regime.

Figure 8: Stability of the numerical solutions in the $(\lambda\rho,\rho)$
plane. Brighter colors denote solutions that tip over more quickly. The black
region shows the parameter range where solutions are considered stable.

Returning to the effective description, we consider the effective (small-
angle) frequency $\omega_{\mathrm{eff}}$, defined in eq. (46), setting $E_{\rm
tot}=0$ in accordance with the initial conditions used to produce Fig. 8. The
lower boundary on $\rho$ for given $\lambda\rho$ corresponds to the stability
condition which was provided by Kapitza at leading power ($\rho^{2}>2$). The
inverted pendulum becomes unstable below that boundary because the driving
amplitude is not sufficient to uphold the pendulum in the inverted position.
Solving eq. (46) for $\omega_{\mathrm{eff}}^{2}=0$ yields a condition on
$\rho$ corresponding to the generalization of eq. (44) to
$\mathcal{O}(\lambda^{26})$. Once these higher-order corrections are included,
the resulting curves approximate the lower boundary between the stable and
unstable region found from the numerical solution to an excellent accuracy, as
illustrated by the blue family of lines in Fig. 9.

Figure 9: Boundary curves obtained from the stability condition
$\omega_{\mathrm{eff}}^{2}>0$ (blue curves) and the convergence condition (48)
(red curves) with increasing accuracy in the $\lambda$ expansion. The shading
of the different curves indicates the value of $n$ used in eq. (48), with the
most saturated color corresponding to the highest $n$. The dark gray area is
the stable region.

The upper boundary of the stability region in Fig. 8 coincides with the value
of $\lambda\rho$ (given $\rho$), for which the $\lambda$ expansion of
$\omega_{\mathrm{eff}}^{2}$ ceases to converge. We can estimate this value by
requiring that the highest available term in the series (46) is smaller than
the previous one, which implies

$\lambda<\sqrt{c_{\omega}^{(2n)}/c_{\omega}^{(2n+2)}}\,,$ (48)

and quantifies the boundary of the region, where $\lambda$ can no longer be
considered small enough to justify the effective treatment. For illustration,
we only show the algebraic result at next-to-leading order ($n=1$),

$\displaystyle\lambda<\frac{4\sqrt{\rho^{2}-2}}{\rho\sqrt{25\rho^{2}-64}}\,,$
(49)

and work with numerical solutions at higher orders. Once these higher-order
results are included, the resulting curves approximate the upper boundary
between the stable and unstable region found from the numerical solution to an
excellent accuracy with increasing approximation order $n$, as illustrated by
the red family of lines in Fig. 9. The least saturated curve refers to $n=1$
and $n=13$ to the most saturated one. Since the effective Lagrangian approach
assumes $\lambda,\lambda\rho\ll 1$, one expects that the boundary lines
obtained from $\omega_{\rm eff}^{2}=0$ (blue) and (48) (red) become unreliable
for too large $\lambda\rho$, hence we do not plot them for $\lambda\rho>1$.

We can therefore conclude that the upright-equilibrium stability boundaries
can be reliably inferred from the condition $\omega_{\rm eff}^{2}=0$ and from
the convergence of the series expansion of $\omega_{\rm eff}^{2}$ in the
effective potential method.

The above findings are an interesting way to study the stability of the
Mathieu functions (i.e. the solutions to the linearized equation of motion of
the system) using effective theories. Another approach which has been adopted
in the literature Blackburn  _et al._ (1992); Mond  _et al._ (1993); Insperger
and Stépán (2002, 2003); Butikov (2011, 2018); Kovacic  _et al._ (2018),
relates the existence of stable, oscillatory solutions to the frequency
separation. In our framework this approach can be stated as follows: When we
assume that the mode functions $\Phi_{n}$ and $\Psi_{n}$ vary with time only
over the time scale of $\omega_{\mathrm{eff}}$, such that the derivatives
scale as powers of $\lambda$, we implicitly make the assumption that
$\omega_{\mathrm{eff}}\sim\lambda\omega=\omega_{0}$. However, as can be seen
from (46), this is no longer the case when $\rho$ becomes increasingly large,
which corresponds to the situation when the amplitude $h$ of the pivot point
oscillation is no longer of order $\lambda$ relative to the length of the
pendulum. From the physical point of view, the criterion for mode separation
is not $\omega\gg\omega_{0}$ but $\omega\gg\omega_{\mathrm{eff}}$. When
$\lambda\rho$ is large enough such that

$\omega_{1}^{-}\equiv\omega-\omega_{\mathrm{eff}}\sim\omega_{\mathrm{eff}}\,,$
(50)

the effective treatment based on fast and slow frequency separation is
certainly no longer applicable. Equating
$\omega_{1}^{-}=\omega_{\mathrm{eff}}$ and inserting the leading-order
expression for $\omega_{\mathrm{eff}}=\lambda\omega\sqrt{\rho^{2}/2-1}$ as an
estimate, yields an approximation for when this occurs, namely at

$\displaystyle\lambda_{\rm
crit}\gtrsim\frac{1}{\sqrt{2}}(\rho^{2}-2)^{-1/2}\,,$ (51)

which at $\rho=2$ amounts to $\lambda_{\mathrm{crit}}=1/2$.

This method provides yet another way of identifying parameter regions in which
the effective description cannot converge. Using our
$\mathcal{O}(\lambda^{28})$-expression for $\omega_{\mathrm{eff}}$, and
solving the relation $\omega_{1}^{-}=\omega_{\mathrm{eff}}$ for $\lambda$
numerically for a given $\rho$, we find that the resulting curves are
virtually indistinguishable from the red lines in Fig. 9 (and therefore not
drawn), although their algebraic expressions are not identical.

##  V Conclusion

In this paper, we applied the effective Lagrangian method beyond the leading
approximation to a classical problem, the driven inverted pendulum with
driving frequency $\omega$ much higher than the one of the unperturbed
pendulum, $\omega_{0}$. Using methods inspired by scale separation and power-
counting in effective quantum field theories, we derived a general algorithm,
valid to any order in the scale separation parameter
$\lambda=\omega_{0}/\omega$, in which the systematic expansion is constructed.
The approach employs a modulated Fourier mode expansion and generalizes the
procedure found in textbooks Landau and Lifshitz (1960), which is based on
averaging high frequency modes.

The effective Lagrangian obtained at higher orders in $\lambda$ is time-
independent, but exhibits interesting properties that would be considered
unconventional and even problematic for a general mechanical system, as it
depends on high powers of the derivative of the generalized coordinate, as
well as on its higher derivatives. These two properties translate to the
system not having a potential of the form $V(\Theta)$, as well as possibly
having an unbounded Hamiltonian. However, the problematic terms appear only in
the Lagrangian at higher order in the expansion parameter. This allows us to
treat them as perturbative corrections and in turn eliminate them in such a
way that the full system is canonical and has an effective potential $V_{\rm
eff}(\Theta,E_{\rm tot})$. The effective potential depends however on the
total energy of the system. By energy conservation following from the time-
independence of the effective Lagrangian, this feature can be rephrased as a
dependence on the initial conditions. The initial conditions themselves need
to be carefully matched between the effective system with generalized
coordinate $\Theta$ and the full system with coordinate $\theta$.

We computed the effective potential explicitly including the first three
corrections to the previously known leading expression, and the curvature of
the potential in the upright equilibrium position to order $\lambda^{28}$,
corresponding to the first 14 terms terms of the series.

By studying the effective potential at such high orders in the expansion
parameters, we estimated the parameter regions of stable trajectories to high
accuracy. This method allows us to determine the parameter-space region in
which parametric resonance occurs, the effective Lagrangian method breaks
down, and the upright equilibrium is unstable. We find that this approach
yields comparable results to the one found in previous studies from imposing
sufficient separation of frequencies.

## Acknowledgments

M.B. thanks A. Penin for making him aware of the inverted pendulum many years
ago. The software FORM Vermaseren (2000) was used for the high-order
calculations. This work has been supported in part by the Excellence Cluster
ORIGINS funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy –EXC-2094 –390783311.

##  Appendix A Higher-Order Effective Potential

We give the subleading expressions for the effective potentials
$V_{\mathrm{eff}}^{(n)}$. For notational convenience, we introduce the
rescaled energy expansion coefficients
$k_{n}=E^{(n)}_{\mathrm{tot}}/\left(\lambda^{n}\omega^{2}\right)$, where
$E^{(n)}_{\mathrm{tot}}$ again refers to the $\lambda^{n}$ contribution of the
total energy. We obtain the following potentials up to the third subleading
correction:

$\displaystyle\frac{V_{\mathrm{eff}}^{(2)}}{\lambda^{2}\omega^{2}}\ =$
$\displaystyle\
c_{\Theta}-1+\rho^{2}\left(-\frac{1}{8}c_{2\Theta}+\frac{1}{8}\right)$
$\displaystyle\frac{V_{\mathrm{eff}}^{(4)}}{\lambda^{4}\omega^{2}}\ =$
$\displaystyle\
\rho^{2}\left(\frac{17}{16}c_{\Theta}-\frac{3}{4}c_{2\Theta}+\frac{7}{16}c_{3\Theta}-\frac{3}{4}+k_{2}\left[-\frac{3}{4}c_{2\Theta}+\frac{3}{4}\right]\right)$
$\displaystyle+\rho^{4}\left(\frac{1}{32}c_{2\Theta}-\frac{29}{512}c_{4\Theta}+\frac{13}{512}\right)$
$\displaystyle\frac{V_{\mathrm{eff}}^{(6)}}{\lambda^{6}\omega^{2}}\ =$
$\displaystyle\
\rho^{2}\left(\frac{87}{4}c_{\Theta}-\frac{25}{2}c_{2\Theta}+\frac{33}{4}c_{3\Theta}-\frac{27}{8}c_{4\Theta}-\frac{113}{8}\right.$
$\displaystyle\qquad\quad+\left.k_{2}\left[\frac{87}{4}c_{\Theta}-5c_{2\Theta}+\frac{33}{4}c_{3\Theta}-25\right]+k_{2}^{2}\left[-\frac{5}{2}c_{2\Theta}+\frac{5}{2}\right]+k_{4}\left[-\frac{3}{4}c_{2\Theta}+\frac{3}{4}\right]\right)$
$\displaystyle+\rho^{4}\left(\frac{2167}{2048}c_{\Theta}-\frac{123}{32}c_{2\Theta}+\frac{10665}{4096}c_{3\Theta}-\frac{629}{256}c_{4\Theta}+\frac{6345}{4096}c_{5\Theta}+\frac{279}{256}\right.$
$\displaystyle\qquad\quad+\left.k_{2}\left[-\frac{87}{32}c_{2\Theta}-\frac{629}{256}c_{4\Theta}+\frac{1325}{256}\right]\right)$
$\displaystyle+\rho^{6}\left(\frac{20183}{73728}c_{2\Theta}+\frac{243}{2048}c_{4\Theta}-\frac{1143}{8192}c_{6\Theta}-\frac{4661}{18432}\right)$
$\displaystyle\frac{V_{\mathrm{eff}}^{(8)}}{\lambda^{8}\omega^{2}}\ =$
$\displaystyle\
\rho^{2}\left(\frac{1981}{4}c_{\Theta}-371c_{2\Theta}+\frac{1989}{8}c_{3\Theta}-\frac{605}{4}c_{4\Theta}+\frac{385}{8}c_{5\Theta}-\frac{1079}{4}\right.$
$\displaystyle\qquad\quad+\left.k_{2}\left[\frac{1211}{2}c_{\Theta}-385c_{2\Theta}+\frac{429}{2}c_{3\Theta}-\frac{605}{4}c_{4\Theta}-\frac{1135}{4}\right]\right.$
$\displaystyle\qquad\quad+\left.k_{2}^{2}\left[\frac{1211}{4}c_{\Theta}-21c_{2\Theta}+\frac{429}{4}c_{3\Theta}-389\right]+k_{2}^{3}\left[-7c_{2\Theta}+7\right]\right.$
$\displaystyle\qquad\quad+\left.k_{4}\left[\frac{87}{4}c_{\Theta}-5c_{2\Theta}+\frac{33}{4}c_{3\Theta}-25\right]+k_{2}k_{4}\left[-5c_{2\Theta}+5\right]+k_{6}\left[-\frac{3}{4}c_{2\Theta}+\frac{3}{4}\right]\right)$
$\displaystyle+\rho^{4}\left(\frac{15903}{128}c_{\Theta}-\frac{7836379}{32768}c_{2\Theta}+\frac{499467}{2048}c_{3\Theta}-\frac{2977261}{16384}c_{4\Theta}\right.$
$\displaystyle\qquad\quad+\left.\frac{261525}{2048}c_{5\Theta}-\frac{1557909}{32768}c_{6\Theta}-\frac{449115}{16384}\right.$
$\displaystyle\qquad\quad+\left.k_{2}\left[\frac{11727}{128}c_{\Theta}-283c_{2\Theta}+\frac{474123}{2048}c_{3\Theta}-\frac{4627}{32}c_{4\Theta}+\frac{261525}{2048}c_{5\Theta}-\frac{2973}{128}\right]\right.$
$\displaystyle\qquad\quad+\left.k_{2}^{2}\left[-\frac{253}{2}c_{2\Theta}-\frac{4627}{64}c_{4\Theta}+\frac{12723}{64}\right]+k_{4}\left[-\frac{87}{32}c_{2\Theta}-\frac{629}{256}c_{4\Theta}+\frac{1325}{256}\right]\right)$
$\displaystyle+\rho^{6}\left(\frac{44553629}{5308416}c_{\Theta}+\frac{2503511}{110592}c_{2\Theta}-\frac{39310133}{5308416}c_{3\Theta}-\frac{81293}{9216}c_{4\Theta}\right.$
$\displaystyle\qquad\quad+\left.\frac{1665323}{196608}c_{5\Theta}-\frac{71277}{4096}c_{6\Theta}+\frac{749151}{65536}c_{7\Theta}-\frac{478423}{27648}\right.$
$\displaystyle\qquad\quad+\left.k_{2}\left[\frac{3386735}{110592}c_{2\Theta}-\frac{47327}{9216}c_{4\Theta}-\frac{71277}{4096}c_{6\Theta}-\frac{223583}{27648}\right]\right)$
$\displaystyle+\rho^{8}\left(-\frac{7968835}{2654208}c_{2\Theta}+\frac{9568931}{10616832}c_{4\Theta}+\frac{333907}{294912}c_{6\Theta}-\frac{1622357}{2097152}c_{8\Theta}+\frac{295983029}{169869312}\right)$

##  Appendix B Higher-Order Initial Conditions

We provide the matched initial condition to the third subleading power. For
simplicity we set the initial time $t_{i}=0$ and define
$\theta_{i}\equiv\theta(0)$. We also define the rescaled initial velocity of
the full system $\nu_{i}\equiv\frac{\dot{\theta}(0)}{\lambda\omega}$. The
expansion of the initial conditions for the generalized coordinate $\Theta(t)$
of the effective description takes the form

$\displaystyle\Theta(0)$
$\displaystyle=\theta_{i}+\sum_{n=1}\Theta_{i}^{(2n)}\,,$
$\displaystyle\dot{\Theta}(0)$
$\displaystyle=\dot{\theta}_{i}-\lambda\rho\omega
s_{\theta_{i}}+\sum_{n=1}\dot{\Theta}_{i}^{(2n+1)}\,,$

where we defined the $\mathcal{O}(\lambda^{2n})$ contribution to the effective
initial condition $\Theta(0)$ as $\Theta_{i}^{(2n)}$, similarly for its
derivative. In the following expressions for the expansion coefficients, we
use the abbreviations $\nu_{i}\equiv\frac{\dot{\theta}(0)}{\lambda\omega}$ and
$\mu_{i}^{(2n)}\equiv\dot{\Theta}_{i}^{(2n+1)}/(\lambda\omega)$:

$\displaystyle\frac{\Theta_{i}^{(2)}}{\lambda^{2}}\ =$ $\displaystyle\ -2\rho
c_{\theta_{i}}\nu_{i}+\frac{17}{16}\rho^{2}s_{2\theta_{i}}$
$\displaystyle\frac{\mu_{i}^{(2)}}{\lambda^{2}}\ =$ $\displaystyle\
\rho\biggl{(}s_{2\theta_{i}}-\nu_{i}^{2}s_{\theta_{i}}\biggr{)}+\rho^{2}\left(\nu_{i}\left[-\frac{3}{8}c_{2\theta_{i}}+\frac{7}{4}\right]\right)+\rho^{3}\left(-\frac{53}{48}s_{\theta_{i}}-\frac{5}{16}s_{3\theta_{i}}\right)$
$\displaystyle\frac{\Theta_{i}^{(4)}}{\lambda^{4}}\ =$ $\displaystyle\
\rho\biggl{(}\nu_{i}\left[11c_{2\theta_{i}}-3\right]-4\nu_{i}^{3}c_{\theta_{i}}\biggr{)}+\rho^{2}\left(\frac{947}{128}s_{\theta_{i}}-\frac{865}{128}s_{3\theta_{i}}+\frac{45}{8}\nu_{i}^{2}s_{2\theta_{i}}\right)$
$\displaystyle+\rho^{3}\left(\nu_{i}\left[-\frac{409}{108}c_{\theta_{i}}-\frac{3}{2}c_{3\theta_{i}}\right]\right)+\rho^{4}\left(-\frac{1529}{3456}s_{2\theta_{i}}+\frac{4039}{2048}s_{4\theta_{i}}\right)$
$\displaystyle\frac{\mu_{i}^{(4)}}{\lambda^{4}}\ =$ $\displaystyle\
\rho\biggl{(}s_{\theta_{i}}-3s_{3\theta_{i}}+11\nu_{i}^{2}s_{2\theta_{i}}-\nu_{i}^{4}s_{\theta_{i}}\biggr{)}$
$\displaystyle+\rho^{2}\left(\nu_{i}\left[-\frac{2067}{128}c_{\theta_{i}}+\frac{899}{128}c_{3\theta_{i}}\right]+\nu_{i}^{3}\left(-\frac{5}{4}c_{2\theta_{i}}+\frac{11}{2}\right)\right)$
$\displaystyle+\rho^{3}\left(\frac{21961}{3456}s_{2\theta_{i}}+\frac{1319}{768}s_{4\theta_{i}}+\nu_{i}^{2}\left[-\frac{5867}{432}s_{\theta_{i}}-\frac{75}{16}s_{3\theta_{i}}\right]\right)$
$\displaystyle+\rho^{4}\left(\nu_{i}\left[\frac{3361}{1728}c_{2\theta_{i}}-\frac{1785}{512}c_{4\theta_{i}}+\frac{42833}{6912}\right]\right)$
$\displaystyle+\rho^{5}\left(-\frac{10507}{34560}s_{\theta_{i}}-\frac{66475}{27648}s_{3\theta_{i}}+\frac{139}{3072}s_{5\theta_{i}}\right)$
$\displaystyle\frac{\Theta_{i}^{(6)}}{\lambda^{6}}\ =$ $\displaystyle\
\rho\biggl{(}\nu_{i}\left[58c_{\theta_{i}}-90c_{3\theta_{i}}\right]+\nu_{i}^{3}\left[114c_{2\theta_{i}}-10\right]-6\nu_{i}^{5}c_{\theta_{i}}\biggr{)}$
$\displaystyle+\rho^{2}\left(-\frac{19437}{256}s_{2\theta_{i}}+\frac{28125}{512}s_{4\theta_{i}}+\nu_{i}^{2}\left[\frac{100257}{512}s_{\theta_{i}}-\frac{86739}{512}s_{3\theta_{i}}\right]+\frac{259}{16}\nu_{i}^{4}s_{2\theta_{i}}\right)$
$\displaystyle+\rho^{3}\left(\nu_{i}\left[\frac{5098439}{31104}c_{2\theta_{i}}-\frac{33565}{13824}c_{4\theta_{i}}-\frac{4253917}{41472}\right]\right.$
$\displaystyle\qquad\quad\left.+\nu_{i}^{3}\left[-\frac{13963}{486}c_{\theta_{i}}-\frac{359}{6}c_{3\theta_{i}}\right]\right)$
$\displaystyle+\rho^{4}\left(\frac{480145957}{7962624}s_{\theta_{i}}+\frac{8066321}{1990656}s_{3\theta_{i}}-\frac{29580613}{884736}s_{5\theta_{i}}\right.$
$\displaystyle\qquad\quad\left.+\nu_{i}^{2}\left[-\frac{212263}{2592}s_{2\theta_{i}}+\frac{53107}{512}s_{4\theta_{i}}\right]\right)$
$\displaystyle+\rho^{5}\left(\nu_{i}\left[\frac{95375231}{1944000}c_{\theta_{i}}-\frac{2177353}{20736}c_{3\theta_{i}}+\frac{125921}{3456}c_{5\theta_{i}}\right]\right)$
$\displaystyle+\rho^{6}\left(-\frac{39063010973}{1327104000}s_{2\theta_{i}}+\frac{10440595}{497664}s_{4\theta_{i}}+\frac{1506521}{3538944}s_{6\theta_{i}}\right)$
$\displaystyle\frac{\mu_{i}^{(6)}}{\lambda^{6}}\ =$ $\displaystyle\
\rho\biggl{(}-18s_{2\theta_{i}}+17s_{4\theta_{i}}+\nu_{i}^{2}\left[29s_{\theta_{i}}-135s_{3\theta_{i}}\right]+57\nu_{i}^{4}s_{2\theta_{i}}-\nu_{i}^{6}s_{\theta_{i}}\biggr{)}$
$\displaystyle+\rho^{2}\left(\nu_{i}\left[\frac{23917}{128}c_{2\theta_{i}}-\frac{13501}{128}c_{4\theta_{i}}-\frac{117}{4}\right]\right.$
$\displaystyle\qquad\quad+\left.\nu_{i}^{3}\left[-\frac{103179}{512}c_{\theta_{i}}+\frac{41811}{512}c_{3\theta_{i}}\right]+\nu_{i}^{5}\left[-\frac{21}{8}c_{2\theta_{i}}+\frac{45}{4}\right]\right)$
$\displaystyle+\rho^{3}\left(\frac{4069973}{62208}s_{\theta_{i}}-\frac{1231805}{20736}s_{3\theta_{i}}-\frac{2369}{432}s_{5\theta_{i}}\right.$
$\displaystyle\qquad\quad+\left.\nu_{i}^{2}\left[\frac{29980835}{124416}s_{2\theta_{i}}+\frac{2973521}{27648}s_{4\theta_{i}}\right]+\nu_{i}^{4}\left[-\frac{113107}{1944}s_{\theta_{i}}-\frac{419}{8}s_{3\theta_{i}}\right]\right)$
$\displaystyle+\rho^{4}\left(\nu_{i}\left[-\frac{725855525}{7962624}c_{\theta_{i}}-\frac{174397183}{1990656}c_{3\theta_{i}}+\frac{98781193}{884736}c_{5\theta_{i}}\right]\right.$
$\displaystyle\qquad\quad+\left.\nu_{i}^{3}\left[\frac{418565}{3888}c_{2\theta_{i}}-\frac{10511}{128}c_{4\theta_{i}}+\frac{1025965}{15552}\right]\right)$
$\displaystyle+\rho^{5}\left(-\frac{60743917147}{1990656000}s_{2\theta_{i}}+\frac{5150996737}{79626240}s_{4\theta_{i}}-\frac{107574157}{8847360}s_{6\theta_{i}}\right.$
$\displaystyle\qquad\quad+\left.\nu_{i}^{2}\left[\frac{485775337}{7776000}s_{\theta_{i}}-\frac{15323549}{82944}s_{3\theta_{i}}+\frac{283495}{27648}s_{5\theta_{i}}\right]\right)$
$\displaystyle+\rho^{6}\left(\nu_{i}\left[\frac{233925934847}{1990656000}c_{2\theta_{i}}-\frac{19863253}{497664}c_{4\theta_{i}}-\frac{32679059}{1769472}c_{6\theta_{i}}-\frac{1259207923}{31104000}\right]\right)$
$\displaystyle+\rho^{7}\left(\frac{183455603621}{5573836800}s_{\theta_{i}}-\frac{67196055007}{3981312000}s_{3\theta_{i}}-\frac{234812363}{31850496}s_{5\theta_{i}}+\frac{63757471}{17694720}s_{7\theta_{i}}\right)$

##  Appendix C Harmonic Frequency Expressions

We give the expansion coefficients for the effective harmonic frequency
$\omega_{\rm eff}^{2}$ defined in (46). We calculated the first 14 terms of
the expansion, which corresponds to $\lambda^{28}$. The coefficients are
polynomials in $\rho$ with rational coefficients.

$\displaystyle c_{\omega}^{(2)}=$ $\displaystyle\frac{1}{2}\rho^{2}-1$
$\displaystyle c_{\omega}^{(4)}=$
$\displaystyle\frac{25}{32}\rho^{4}-2\rho^{2}$ $\displaystyle
c_{\omega}^{(6)}=$
$\displaystyle\frac{1169}{576}\rho^{6}-\frac{273}{32}\rho^{4}+8\rho^{2}$
$\displaystyle c_{\omega}^{(8)}=$
$\displaystyle\frac{16824665}{2654208}\rho^{8}-\frac{94831}{2592}\rho^{6}+\frac{2049}{32}\rho^{4}-32\rho^{2}$
$\displaystyle c_{\omega}^{(10)}=$
$\displaystyle\frac{104550461873}{4777574400}\rho^{10}-\frac{15301534105}{95551488}\rho^{8}+\frac{18943067}{46656}\rho^{6}-\frac{13057}{32}\rho^{4}+128\rho^{2}$
$\displaystyle c_{\omega}^{(12)}=$
$\displaystyle\frac{1383860829361699}{17199267840000}\rho^{12}-\frac{510599414365027}{716636160000}\rho^{10}+\frac{2704847264795}{1146617856}\rho^{8}-\frac{187198453}{52488}\rho^{6}+\frac{75777}{32}\rho^{4}-512\rho^{2}$
$\displaystyle c_{\omega}^{(14)}=$
$\displaystyle\frac{29289023958538918009}{94810963968000000}\rho^{14}-\frac{12419023504555535053}{3869835264000000}\rho^{12}+\frac{16897194776220398401}{1289945088000000}\rho^{10}-\frac{3275894575046105}{123834728448}\rho^{8}$
$\displaystyle+\frac{102745082885}{3779136}\rho^{6}-\frac{413697}{32}\rho^{4}+2048\rho^{2}$
$\displaystyle c_{\omega}^{(16)}=$
$\displaystyle\frac{83857909128685125912504026393}{68504182964040499200000000}\rho^{16}-\frac{488045135705765455513091071}{33449308087910400000000}\rho^{14}+\frac{985874975154709612644967}{13931406950400000000}\rho^{12}$
$\displaystyle-\frac{155733295202139312167933}{870712934400000000}\rho^{10}+\frac{1114567661360939425}{4458050224128}\rho^{8}-\frac{3200742128665}{17006112}\rho^{6}+\frac{2162689}{32}\rho^{4}-8192\rho^{2}$
$\displaystyle c_{\omega}^{(18)}=$
$\displaystyle\frac{120089176174908625267123147897106657}{24168275749713488117760000000000}\rho^{18}-\frac{269168508120402223082153608059967051}{4028045958285581352960000000000}\rho^{16}$
$\displaystyle+\frac{92279593282170513811138202650217}{245852414446141440000000000}\rho^{14}-\frac{3571425412964887484418187991}{3134566563840000000000}\rho^{12}+\frac{2108768571265773123745227511}{1044855521280000000000}\rho^{10}$
$\displaystyle-\frac{112207057270371335395}{53496602689536}\rho^{8}+\frac{371789779789295}{306110016}\rho^{6}-\frac{10944513}{32}\rho^{4}+32768\rho^{2}$
$\displaystyle c_{\omega}^{(20)}=$
$\displaystyle\frac{175263403904227897478594054729227794909043}{8526567684498918607945728000000000000}\rho^{20}-\frac{131243335963394451160758386540543920310989}{426328384224945930397286400000000000}\rho^{18}$
$\displaystyle+\frac{279232015345030240850013025618323421714447}{142109461408315310132428800000000000}\rho^{16}-\frac{45296196240036131201719335763017237187}{6505254886244902502400000000000}\rho^{14}$
$\displaystyle+\frac{33739061859510881881535523290359}{2256887925964800000000000}\rho^{12}-\frac{9355356107413400058873181771417}{470184984576000000000000}\rho^{10}+\frac{93058155781360081771505}{5777633090469888}\rho^{8}$
$\displaystyle-\frac{5116060734204095}{688747536}\rho^{6}+\frac{54001665}{32}\rho^{4}-131072\rho^{2}$
$\displaystyle c_{\omega}^{(22)}=$
$\displaystyle\frac{282847525631781705724885459394868767290384277872373}{3275900483130336930037862341017600000000000000}\rho^{22}-\frac{19290040218535541559219833629556583270732139031503}{13536778855910483181974637772800000000000000}\rho^{20}$
$\displaystyle+\frac{13788286801022032955390760875602758321381259880239}{1353677885591048318197463777280000000000000}\rho^{18}-\frac{3103818637959156164712319406632069993514674355453}{75204326977280462122081320960000000000000}\rho^{16}$
$\displaystyle+\frac{7470836385204179119168031963532059214635417}{71720435120850050088960000000000000}\rho^{14}-\frac{430039973727914354734435397867644739}{2538998916710400000000000000}\rho^{12}$
$\displaystyle+\frac{449916002199871703719766958103238663}{2538998916710400000000000000}\rho^{10}-\frac{24019279122807117265238905}{207994791256915968}\rho^{8}$
$\displaystyle+\frac{1079536847309190425}{24794911296}\rho^{6}-\frac{261095425}{32}\rho^{4}+524288\rho^{2}$
$\displaystyle c_{\omega}^{(24)}=$
$\displaystyle\frac{462268189407622980070187142794131930972283406622989689224253}{1258598344898288944251002681129089105920000000000000000}\rho^{24}$
$\displaystyle-\frac{86826106905618634905556618093828366949132731009071000532877}{13110399426023843169281277928428011520000000000000000}\rho^{22}$
$\displaystyle+\frac{93847416638175046132782637159577187302841621594668517389}{1790915842636956924975244577341440000000000000000}\rho^{20}$
$\displaystyle-\frac{128387128420533705367916342309325844669854222563018198989}{537274752791087077492573373202432000000000000000}\rho^{18}$
$\displaystyle+\frac{1022077886396251550445416469859905138490704725080596389}{1474004808754697057592793890816000000000000000}\rho^{16}$
$\displaystyle-\frac{6248781427949306780090443236918867305119283092831}{4685735094562203272478720000000000000000}\rho^{14}+\frac{5229382962963418635872640393212715850861}{3046798700052480000000000000000}\rho^{12}$
$\displaystyle-\frac{46248045973513272796578928600044367951}{31737486458880000000000000000}\rho^{10}+\frac{217332854394611010157695955}{277326388342554624}\rho^{8}-\frac{27507355568072365315}{111577100832}\rho^{6}$
$\displaystyle+\frac{1241513985}{32}\rho^{4}-2097152\rho^{2}$ $\displaystyle
c_{\omega}^{(26)}=$
$\displaystyle\frac{129038512518653608908013841918133236441640681858146168536514742761447177}{81720368652080691243763491149613041176717295616000000000000000000}\rho^{26}$
$\displaystyle-\frac{7466672475271258206411032946255392638581353745701451188541219236031989}{241776238615623346875039914643825565611589632000000000000000000}\rho^{24}$
$\displaystyle+\frac{5398845519120437529448129129031926002295887567322877054887392159172379}{20148019884635278906253326220318797134299136000000000000000000}\rho^{22}$
$\displaystyle-\frac{23168943557410693360731764696486962814988326646666342872117918209}{17059547950622596884544189745923620864000000000000000000}\rho^{20}$
$\displaystyle+\frac{15199963329184432993009020771898106646677055890472950214498023017}{3411909590124519376908837949184724172800000000000000000}\rho^{18}$
$\displaystyle-\frac{12869476934636564747653863325226820557044411384935409167629513}{1300072241321642804796844211699712000000000000000000}\rho^{16}$
$\displaystyle+\frac{14089529769328588448840379000605568208323725234457151441}{929884129515869239423401984000000000000000000}\rho^{14}$
$\displaystyle-\frac{10915670389386472594378595790158606675747899}{685529707511808000000000000000000}\rho^{12}+\frac{2567346427624746880813248711133145887270177}{228509902503936000000000000000000}\rho^{10}$
$\displaystyle-\frac{1370308410719375473886395601225}{269561249468963094528}\rho^{8}+\frac{2724400629069976643459}{2008387814976}\rho^{6}-\frac{5821693953}{32}\rho^{4}+8388608\rho^{2}$
$\displaystyle c_{\omega}^{(28)}=$
$\displaystyle\frac{36345946021027445341629359633703104737621781477317320811023948708058852109306767193}{5306076143912027105009858570822141413312851568709389516800000000000000000000}\rho^{28}$
$\displaystyle-\frac{27369000728578716214966200477234554367256981455962616682322216498296622080436149283}{189502719425429539464637806100790764761173270311049625600000000000000000000}\rho^{26}$
$\displaystyle+\frac{9046641423542327760184622591380741775142835883666872154875740353786511874617599}{6635016961080828383622345369587576231965731953049600000000000000000000}\rho^{24}$
$\displaystyle-\frac{6306250910017840861488192420418119232469566746100719153263715392889596030343689}{829377120135103547952793171198447028995716494131200000000000000000000}\rho^{22}$
$\displaystyle+\frac{215301979372622936192888758236579557739971628689805217012783163354456239}{7738210950402409946829244468750954423910400000000000000000000}\rho^{20}$
$\displaystyle-\frac{135744408063298052316132454385830653034704388562887895544837073042367557}{1934552737600602486707311117187738605977600000000000000000000}\rho^{18}$
$\displaystyle+\frac{12250032304270105368254748800235591753615443417621184518084360837597}{98285461443916196042641422404498227200000000000000000000}\rho^{16}$
$\displaystyle-\frac{3209650102469216338239811956334222922347371427779307065627301}{20503945055824916729286013747200000000000000000000}\rho^{14}$
$\displaystyle+\frac{339268309048445977552318981024838353042166016577}{2467906947042508800000000000000000000}\rho^{12}-\frac{76004060694084228949294349897409420899453027821}{925465105140940800000000000000000000}\rho^{10}$
$\displaystyle+\frac{308135893627892285068404828937105}{9704204980882671403008}\rho^{8}-\frac{8230657846197335693237}{1129718145924}\rho^{6}+\frac{26977763329}{32}\rho^{4}-33554432\rho^{2}$

## References

  * Stephenson (1908) A. Stephenson, Lond. Edinb. Dublin Philos. Mag. J. Sci. 15, 233 (1908).
  * Kapitza (1951) P. Kapitza, J. Exp. Theor. Phys. 21, 588 (1951).
  * Landau and Lifshitz (1960) L. D. Landau and Lifshitz, in  _Mechanics_ , Vol. 1 (Pergamon Press, Oxford, 1960) Chap. 30.
  * Butikov (2001) E. Butikov, Am. J. Phys. 69 (2001).
  * Coleman and Weinberg (1973) S. R. Coleman and E. J. Weinberg, Phys. Rev. D7, 1888 (1973).
  * Buttazzo  _et al._ (2013) D. Buttazzo, G. Degrassi, P. P. Giardino, G. F. Giudice, F. Sala, A. Salvio, and A. Strumia, JHEP 12, 089 (2013).
  * Rahav  _et al._ (2003a) S. Rahav _et al._ , Phys. Rev. Lett. 91, 110404 (2003a).
  * Grundy (2019) R. E. Grundy, Q. J. Mech. Appl. Math. 72, 261 (2019).
  * Maggia  _et al._ (2020) M. Maggia _et al._ , Nonlinear Dyn. 99, 813 (2020).
  * Kuehn (2015) C. Kuehn,  _Multiple Time Scale Dynamics_ (Springer, Berlin, 2015) pp. 277–278.
  * Sanz-Serna (2009) J. M. Sanz-Serna, IMA J. Numer. Anal. 29, 595 (2009).
  * Butikov (2011) E. I. Butikov, J. Phys. A Math. Theor. 44, 295202 (2011).
  * Butikov (2018) E. I. Butikov, Am. J. Phys. 86, 257 (2018).
  * Blackburn  _et al._ (1992) J. A. Blackburn, H. Smith, and N. Grønbech-Jensen, Am. J. Phys. 60, 903 (1992).
  * Mond  _et al._ (1993) M. Mond, G. Cederbaum, P. Khan, and Y. Zarmi, J. Sound Vib. 167, 77 (1993).
  * Insperger and Stépán (2002) T. Insperger and G. Stépán, Proc. R. Soc. A: Math. Phys. Eng. 458, 1989 (2002).
  * Insperger and Stépán (2003) T. Insperger and G. Stépán, J. Dyn. Sys., Meas., Control 125, 166 (2003).
  * Jackson (2021) J. D. Jackson,  _Classical electrodynamics_ (John Wiley & Sons, Hoboken, 2021).
  * Weinberg (1979) S. Weinberg, Physica A 96, 327 (1979).
  * Polchinski (1992) J. Polchinski, in  _Theoretical Advanced Study Institute (TASI 92): From Black Holes and Strings to Particles_ (1992) pp. 235–276.
  * Manohar (2020) A. V. Manohar, in  _Les Houches Summer School: EFT in Particle Physics and Cosmology, July 3-28, 2017, published in Les Houches Lect. Notes 108_ (2020).
  * Rahav  _et al._ (2003b) S. Rahav, I. Gilary, and S. Fishman, Phys. Rev. A 68, 013820 (2003b).
  * Ostrogradsky (1850) M. Ostrogradsky, Mem. Acad. St. Petersbourg 6, 385 (1850).
  * Anderson (1973) D. Anderson, J. Phys., A Math. Nucl. Gen. 6, 299 (1973).
  * Cohen and Green (2020) T. Cohen and D. Green, JHEP 12, 041 (2020).
  * Cohen  _et al._ (2021) T. Cohen, D. Green, A. Premkumar, and A. Ridgway, JHEP 09, 159 (2021).
  * Bender and Wu (1969) C. M. Bender and T. T. Wu, Phys. Rev. 184, 1231 (1969).
  * Dyson (1952) F. J. Dyson, Phys. Rev. 85, 631 (1952).
  * Kovacic  _et al._ (2018) I. Kovacic, R. Rand, and S. Mohamed Sah, Appl. Mech. Rev. 70, 020802 (2018).
  * Vermaseren (2000) J. A. Vermaseren, arXiv preprint math-ph/0010025 (2000).