# When action is not least for systems with action-dependent Lagrangians

Joseph Ryan<EMAIL_ADDRESS>Department of Physics, Southern Methodist
University, Dallas, TX 75275, USA

###### Abstract

The dynamics of some non-conservative and dissipative systems can be derived
by calculating the first variation of an action-dependent action, according to
the variational principle of Herglotz. This is directly analogous to the
variational principle of Hamilton commonly used to derive the dynamics of
conservative systems. In a similar fashion, just as the second variation of a
conservative system’s action can be used to infer whether that system’s
possible trajectories are dynamically stable, so too can the second variation
of the action-dependent action be used to infer whether the possible
trajectories of non-conservative and dissipative systems are dynamically
stable. In this paper I show, generalizing earlier analyses of the second
variation of the action for conservative systems, how to calculate the second
variation of the action-dependent action and how to apply it to two physically
important systems: a time-independent harmonic oscillator and a time-dependent
harmonic oscillator.

##  I Introduction

The principle of stationary action is a powerful unifying principle for
physics. The dynamics of any conservative system, whose configuration at a
given time $t$ is described by the generalized coordinates
$\vec{x}(t)=\left(x_{1}(t),x_{2}(t),...,x_{n}(t)\right)^{\rm T}$ and the
generalized velocity $\vec{v}(t)=\frac{d\vec{x}(t)}{dt}$, can be derived by
demanding that the functional

$S=\int_{t_{i}}^{t_{f}}L\left(\vec{x}(t),\vec{v}(t)\right)dt,$ (1)

be stationary with respect to variations of the form

$\vec{x}\left(t,\epsilon\right)=\vec{x}\left(t,0\right)+\epsilon\vec{y}\left(t\right),$
(2)

where $y(t)$ satisfies $y\left(t_{i}\right)=y\left(t_{f}\right)=0$ but is
otherwise arbitrary, and $\epsilon$ is an arbitrary infinitesimal
parameter.Goldstein, Poole, and Safko (2002); Fowles, Grant R. and Cassiday,
George L. (2005) By “stationary”, I mean that

$\frac{dS}{d\epsilon}\bigg{|}_{\epsilon\rightarrow 0}=0.$ (3)

The system will follow the trajectories that make the action stationary as
defined by eq. (3). For this condition to be satisfied, the Lagrangian
$L\left(\vec{x}(t),\vec{v}(t)\right)$ must obey the Euler-Lagrange equations

$\frac{d}{dt}\frac{\partial L}{\partial\dot{x}_{i}}=\frac{\partial L}{\partial
x_{i}},$ (4)

where $x_{i}$ refers to the $i^{\rm th}$ component of the vector $\vec{x}(t)$
and $\dot{x}_{i}:=v_{i}$. The trajectories that make the action stationary
must therefore be the ones that satisfy eq. (4). In non-relativistic
mechanics, the Lagrangian takes the form $L=K-U$, where $K$ is the kinetic
energy of the system and $U$ is the potential energy of the system. Finally,
it is important to note that the Euler-Lagrange equations are valied for any
time $t_{i}\leq t\leq t_{f}$, and not just at the endpoints $t_{i}$ and
$t_{f}$ (this will be important in Sec. II.2).

This variational, energy-based approach to mechanics can be a very powerful
alternative to a force-based approach, and it is more easily generalized to
the quantum domain. Fowles, Grant R. and Cassiday, George L. (2005) As
formulated above, however, it cannot be applied to non-conservative systems.
There has been much discussion of this point in the literature (see, for
example, Refs. 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15) and many
techniques have been proposed for handling non-conservative systems within a
variational framework. One technique is to use a Lagrangian that depends
explicitly on time such that $L=L\left(t,\vec{x},\vec{v}\right)$, to
incorporate the effects of, for example, time-dependent external fields or
time-dependent constraints.Gray, Karl, and Novikov (2004) Dissipation that is
linear in the velocity can be accounted for by using a Lagrangian of the form
$L=e^{kt}\left(K-U\right)$, where $k$ is a constant, although there is a long-
running debate in the literature on whether this produces physically
meaningful energy and momentum functions.111See, for example, Refs. 3; 10; 59;
7; 4; 60; 61; 6; 62; 63; 64 and references therein. Another way to handle non-
conservative and dissipative systems, which does produce physically meaningful
energy and momentum functions, is to use a Lagrangian that is an explicit
function of the action $S$ (it may also be an explicit function of time).Lazo
_et al._ (2018) This approach was developed by Gustav Herglotz,Guenther, R.
B., Guenther, C. M., and Gottsch, J. A. (1996) and it has recently been
explored by several authors (see, for example, Refs. 20; 21; 18; 22; 23; 24;
25; 26; 27). A Lagrangian that depends explicity on $S$ can be used to derive
equations of motion that describe certain types of non-conservative and
dissipative systems, and it contains the standard conservative Lagrangian as a
special case. Some investigators have produced generalized versions of
Noether’s theorems stemming from the symmetries of the Herglotz
Lagrangian,Georgieva, Guenther, and Bodurov (2003); Georgieva and Guenther
(2005); Zhang and Tian (2019) and have used these to discover conservation
laws for non-conservative systems Donchev (2014) (see also Ref. 25 for an
analysis of adiabatic invariants for non-conservative systems, and Ref. 28 for
a geometric approach to the Noether symmetry analysis). Others have broadened
the applicability of the Herglotz Lagrangian to field theory, in both its non-
covariant Georgieva, Guenther, and Bodurov (2003) and its covariantLazo  _et
al._ (2017, 2018) forms, leading to alternative theories of gravity and
cosmology Lazo  _et al._ (2017); Paiva, Lazo, and Zanchin (2021); Sloan
(2021); Caramês  _et al._ (2018)

The stationary action principle, whether it is applied to conservative or to
non-conservative systems, is sometimes called the “least action principle”,
but this can be misleading in some situations. It is not always the case that
the action of a conservative system is minimized, because eq. (3) can pick out
a saddle point as well as a minimum (and, as shown in Ref. 30, the action
never attains a true maximum). To answer the question of whether the system’s
action has a true minimum or a saddle point requires the calculation of the
second variation of the system’s action along a stationary path, and it turns
out to be both interesting and useful to do this. In classical mechanics, the
second variation of a system’s action has a close connection to its dynamical
stability, which I will explain in more detail below.Papastavridis (1982,
1983, 1985) Additionally, in quantum mechanics, the phase of the semiclassical
propagator $K(t_{i},t_{f})$ is determined by the sign of the second variation
of the classical action.222See Ref. 45. Technically the phase is equal to
$-\nu\pi/2$, where $\nu$ is the number of negative eigenvalues of the operator
$\Lambda_{ij}$, which I define in eq. (60). Analyses of the second variation
of the action for conservative, classical systems have been carried out in
Refs. 30; 35; 36; 32, and some studies of the second variation of the action
for non-conservative systems exist (see, for example, Refs. 31; 33; 37; 30)
but as far as I know the dynamical stability properties associated with the
second variation of the Herglotz action have not been investigated.333An
explicit calculation of the second variation of the Herglotz action, similar
to what I derive in Sec. II.2, can be found in Ref. 65. As far as I can tell,
however, the authors of that paper appear to have been concerned only with
proving that the solution of the Herglotz variational problem minimizes the
Herglotz action, and not with investigating the sign change of the second
variation. Additionally, these authors do not appear to make the assumption,
as I do, that the second variation of the Herglotz action vanishes at the
initial time $t_{i}$. My goal, therefore, is to generalize the earlier
analyses of the second variation of the action for conservative systems,
focusing mainly on Refs. 30; 36. I will do this by:

1.) showing that the sign of the second variation of the Herglotz action can
be related to the dynamical stability of non-conservative systems, following
Ref. 30, and

2.) showing that, with a few changes to the relevant equations, the sign of
the second variation of the action can be calculated by way of a linear
eigenvalue equation, in the same manner as for a conservative system
(following Ref. 36).

In the next section I summarize the Herglotz variational principle, starting
with the first variation in Sec. II.1, moving on to the second variation in
Sec. II.2. I then derive the connection between the sign of the second
variation of the Herglotz action and dynamical stability in Sec. III. Sec. IV
describes a method for determining the sign of the second variation of the
Herglotz action, followed by two simple examples in Secs. IV.1 and IV.2.
Throughout this paper I use the Einstein summation convention for repeated
indices. For readers who may wish to skip the lengthy derivations, the
following list provides a summary of the key material:

  * •

Euler-Lagrange equation for a system described by the Herglotz action, eq.
(17).

  * •

Second variation of the Herglotz action, eq. (32), followed by comparison to
the second variation of the standard conservative action in eq. (33).

  * •

Linear eigenvalue equation for the action eigenvalue spectrum of a one-
dimensional system with time-dependent dissipation, eq. (73).

  * •

Action eigenvalue spectrum of a simple harmonic oscillator with time-
independent dissipation, eq. (76).

  * •

Numerical estimate of the lowest order action eigenvalue of a simple harmonic
oscillator with time-dependent dissipation and time-dependent frequency, as a
function of final time $t_{f}$, Fig. 2. Validation of this result, using a
numerical solution of the equations of motion for the harmonic oscillator with
time-dependent dissipation and time-dependent frequency, Fig. 3.

##  II Herglotz variational principle

###  II.1 First variation

To obtain the equation(s) of motion, it is necessary to calculate the
functional derivative, not of an integral equation, but of a differential
equation

$\frac{dS}{dt}=L\left(t,x_{i},\dot{x}_{i},S\right),$ (5)

with the initial condition on the action being
$S(t_{i}=0)=:S_{t_{i}}$.Guenther, R. B., Guenther, C. M., and Gottsch, J. A.
(1996) The reader should note here that the action $S$ defined in eq. (5)
depends on the time $t$, unlike the action $S$ defined in eq. 1, which does
not (assuming $t_{f}$ and $t_{i}$ are not allowed to vary). In the rest of
this paper, $t_{i}=0$. We want $S$ to be stationary with respect to variations
of the form shown in eq. (2), and by “stationary” I mean that

$\frac{dS}{d\epsilon}\bigg{|}_{\epsilon\rightarrow 0,t=t_{f}}=0,$ (6)

which is similar to, but not identical with, the condition that conservative
systems must satisfy. In contrast to the standard conservative action, the
first variation of the action vanishes only at the final time
$t_{f}$.Guenther, R. B., Guenther, C. M., and Gottsch, J. A. (1996) Taking the
functional derivative of eq. (5), we find

$\displaystyle\frac{d}{d\epsilon}\left[\frac{dS}{dt}\right]\bigg{|}_{\epsilon\rightarrow
0}$
$\displaystyle=\frac{d}{dt}\frac{dS}{d\epsilon}\bigg{|}_{\epsilon\rightarrow
0}$ (7) $\displaystyle=\frac{\partial L}{\partial x_{i}}y_{i}+\frac{\partial
L}{\partial\dot{x}_{i}}\dot{y}_{i}+\frac{\partial L}{\partial
S}\frac{dS}{d\epsilon}\bigg{|}_{\epsilon\rightarrow 0},$

where the partial derivatives of $L$ are evaluated in the limit
$\epsilon\rightarrow 0$ (that is, on the trajectory that solves the equations
of motion). In the first line of eq. (7), the operators
$\frac{d}{dt}\frac{d}{d\epsilon}$ commute because $\epsilon$ is independent of
time. If we define

$f:=\frac{dS}{d\epsilon}\bigg{|}_{\epsilon\rightarrow 0},$ (8)

and

$B:=\frac{\partial L}{\partial S},$ (9)

then eq. (7) becomes

$\frac{df}{dt}=\frac{\partial L}{\partial x_{i}}y_{i}+\frac{\partial
L}{\partial\dot{x}_{i}}\dot{y}_{i}+Bf.$ (10)

This equation has the form $\frac{df}{dt}=\left(...\right)+Bf$, suggesting
that the solution should be something like

$f=Ae^{\int B\left(\tau\right)d\tau}+{\rm const},$ (11)

where $A$ is an as yet undetermined function of time. The derivative of $f$ is

$\frac{df}{dt}=\frac{dA}{dt}e^{\int B\left(\tau\right)d\tau}+Bf,$ (12)

and by comparing this equation to eq. (10), we can see that eq. (11) solves
eq. (10) if

$\frac{dA}{dt}:=\left(\frac{\partial L}{\partial x_{i}}y_{i}+\frac{\partial
L}{\partial\dot{x}_{i}}\dot{y}_{i}\right)e^{-\int B(\tau)d\tau}$ (13)

More formally, the solution of eq. (10) is given by

$f(t_{f})=Ae^{\int_{0}^{t_{f}}B(\tau)d\tau}+f_{i},$ (14)

where $f_{i}:=\frac{d}{d\epsilon}S_{t_{i}}\big{|}_{\epsilon\rightarrow 0}$.

Because $S_{t_{i}}$ is independent of $x$, it follows that
$\frac{d^{n}}{d\epsilon^{n}}S_{t_{i}}\big{|}_{\epsilon\rightarrow 0}=0$, where
$n$ is an integer, so $f_{i}=0$.Vermeeren, Bravetti, and Seri (2019) At time
$t_{f}$, we demand $\frac{dS}{d\epsilon}\big{|}_{\epsilon\rightarrow 0}=0$, so
$f(t_{f})$ vanishes and eq. 14 becomes

$0=\int_{0}^{t_{f}}\left(\frac{\partial L}{\partial x_{i}}y_{i}+\frac{\partial
L}{\partial\dot{x}_{i}}\dot{y}_{i}\right)e^{-\int_{0}^{\tau_{2}}B(\tau_{1})d\tau_{1}}d\tau_{2}.$
(15)

Integrating the second term in parentheses in eq. (15) by parts, and using the
fact that $y_{i}(0)=y_{i}(t_{f})=0$ to discard the boundary terms, we find

$\displaystyle 0=$ $\displaystyle\int^{t_{f}}_{0}\frac{\partial L}{\partial
x_{i}}y_{i}e^{-\int_{0}^{\tau_{2}}B(\tau_{1})d\tau_{1}}d\tau_{2}$ (16)
$\displaystyle-\int^{t_{f}}_{0}\left[\frac{d}{dt}\frac{\partial
L}{\partial\dot{x}_{i}}y_{i}e^{-\int_{0}^{\tau_{2}}B(\tau_{1})d\tau_{1}}-B\frac{\partial
L}{\partial\dot{x}_{i}}y_{i}e^{-\int_{0}^{\tau_{2}}B(\tau_{1})d\tau_{1}}\right]d\tau_{2},$

which requires

$\frac{d}{dt}\frac{\partial L}{\partial\dot{x}_{i}}=\frac{\partial L}{\partial
x_{i}}+\frac{\partial L}{\partial S}\frac{\partial L}{\partial\dot{x}_{i}}.$
(17)

Eq. (17) is the generalized Euler-Lagrange equation for a system with non-
conservative dynamics, described by an action-dependent Lagrangian.
Applications of this generalized Euler-Lagrange equation can be found in, for
example, Refs. 18; 21; 22; 26; 23; 20; 40; 29.

Figure 1: Solid and dashed curves of the same color have slightly different
amplitudes (see Table 1), and their points of intersection (kinetic foci) lie
on the dotted black curve (the caustic). Displacements recorded on the
vertical axis are measured in arbitrary units, and the times recorded on the
horizontal axis have been scaled by a factor of $\sqrt{\frac{C}{m}}$.

For conservative systems, the first variation of the action is typically used
to derive the equation(s) of motion, and the second variation of the action
can be used to infer whether the system is dynamically stable, as I mentioned
above.Gutzwiller (1990); Gray and Taylor (2007); Papastavridis (1983) To see
how this works, consider a trajectory $\vec{x}_{0}(t)$ that satisfies a given
system’s equation(s) of motion. As shown in Ref. 30, the second variation of
the action for a conservative system is minimized, for all variations
$\\{\epsilon\vec{y}(t)\\}$ away from the stationary path, as long as
$t_{f}<t_{\rm KF1}$, where $t_{f}$ is the final time and “KF1" stands for
“first kinetic focus". As defined in Refs. 30 and 35, a kinetic focus (also
known as a “conjugate point”Gray and Poisson (2011)) is a point along
$\vec{x}_{0}(t)$ at which a nearby trajectory coalesces with it in the limit
that the difference of the initial velocities of both trajectories vanishes.
More precisely, a kinetic focus is a point in time on $\vec{x}_{0}(t)$ that
satisfies

${\rm det}\left[\frac{\partial\vec{x}_{0}(t)}{\partial\vec{v}_{0}}\right]=0,$
(18)

where $\vec{v}_{0}$ is the initial velocity of the particle following
$\vec{x}_{0}(t)$ (note that, for one-dimensional trajectories, eq. (18)
reduces to $\frac{\partial x_{0}(t)}{\partial v_{0}}=0$). From these
equations, we can see that kinetic foci play an important role in
characterizing the dynamical stability of a given system. In a dynamically
unstable system, two nearby trajectories with similar initial velocities will
move farther apart as time goes on. The possible trajectories spread, or “fan
out” Gutzwiller (1990); Gray and Taylor (2007) in a dynamically unstable
system, whereas the possible trajectories of a system that is dynamically
stable will converge, in the sense described above.Gutzwiller (1990); Gray and
Taylor (2007) When $\vec{x}_{0}(t)$ reaches its first kinetic focus, the
second variation of the system’s action will be minimized for all but one of
the variations away from $\vec{x}_{0}(t)$. This special variation (call it
$\epsilon\vec{y}_{*}(t)$) will be another true trajectory of the system,
because it satisfies the system’s equation(s) of motion.Gray and Taylor
(2007); Gray and Poisson (2011) Further, when $t_{f}>t_{\rm KF1}$, the second
variation of the system’s action will be maximized with respect to
$\epsilon\vec{y}_{*}(t)$, and minimized with respect to all other variations
(until the system reaches another kinetic focus; for more details see Refs.
30; 35). In general, then, the action of a conservative system is minimized if
$t_{f}<t_{\rm KF1}$, and is a saddle point if $t_{f}>t_{\rm KF1}$.Gray and
Taylor (2007); Papastavridis (1983)

Roughly speaking, then, a kinetic focus is a kind of “spacetime focal point"
for trajectories having similar initial velocities. A visual example of
kinetic foci is provided by Fig. 1. In this figure, following the example of
Fig. 4 in Ref. 30, I show several possible trajectories of a quartic
oscillator (that is, an oscillator with $U=\frac{1}{4}Cx^{4}$). As in Refs. 30
and 42, I use approximate solutions of the form

$x\approx A{\rm sin}\left(\omega_{0}t\right),$ (19)

where

$\omega_{0}\approx\sqrt{\frac{3C}{4m}}A.$ (20)

Redefining the time variable to have the form

$\tau:=\sqrt{\frac{C}{m}}t,$ (21)

gives

$x\approx A{\rm sin}\left(\sqrt{\frac{3}{4}}A\tau\right).$ (22)

The trajectories shown in Fig. 1 differ in their amplitudes according to Table
1. The important thing to note here is that solid and dashed curves of the
same color have amplitudes that are only slightly different. Because
$v_{0}=\omega_{0}A$, a slight change in the amplitude of a given trajectory
corresponds to a slight change in that trajectory’s initial velocity. From
Fig. 1, we can see that trajectories having similar initial velocities
eventually intersect; these intersection events are called kinetic foci.444The
black dotted curve in the figure is known as the caustic, and it is an
envelope that divides spacetime into regions that are accessible from the
given initial conditions, and regions that aren’t (see Ref. 30). Generally
speaking, the caustic is composed of the kinetic foci associated with the
family of trajectories that are dynamically allowed within a given potential
(i.e. the family of solutions to the system’s equation or equations of
motion). I will not be concerned with caustics in the rest of this paper,
because the examples I consider in Secs. IV.1 and IV.2, namely harmonic
oscillators with time-independent and time-dependent damping and frequency, do
not exhibit caustics. Harmonic oscillators are special in that trajectories
with arbitrarily large differences in their initial velocities reach the same
kinetic focus, so for these systems the caustic curve collapses to a point
Gray and Taylor (2007). The reader should note that I have only shown the
first kinetic foci associated with the trajectories in Fig. 1, and it is at
these points that the second variation of the action, evaluated on a
particular trajectory, vanishes for the first time.

Quartic oscillator trajectories |  $A$ (arb. units)
---|---
$x_{0}$ | 1
$x_{1}$ | 1.05
$x_{2}$ | 1.5
$x_{3}$ | 1.55
$x_{4}$ | 2
$x_{5}$ | 2.05
Table 1: Trajectories used to generate Fig. 1. The left column lists the
trajectories with their corresponding label subscripts, and the right column
lists the initial values of their amplitudes in arbitrary units.

I will derive a general expression for the second variation of the Herglotz
action in the next section, which I will use, in Sec. III, to generalize the
analysis of Ref. 30.

###  II.2 Second variation

To calculate the second functional derivative of the Herglotz action, we take
the functional derivative of eq. (7):

$\displaystyle\frac{d}{d\epsilon}\left[\frac{d}{d\epsilon}\frac{dS}{dt}\right]\bigg{|}_{\epsilon\rightarrow
0}$
$\displaystyle=\frac{d}{d\epsilon}\left[\frac{d}{dt}\frac{dS}{d\epsilon}\right]\bigg{|}_{\epsilon\rightarrow
0}$ (23)
$\displaystyle=\frac{d}{d\epsilon}\left[\frac{df}{dt}\right]\bigg{|}_{\epsilon\rightarrow
0}$
$\displaystyle=\left[\frac{d\tilde{A}}{d\epsilon}\right]\bigg{|}_{\epsilon\rightarrow
0}+\left[\frac{dB}{d\epsilon}\right]\bigg{|}_{\epsilon\rightarrow
0}f+B\left[\frac{df}{d\epsilon}\right]\bigg{|}_{\epsilon\rightarrow 0},$

where $B$ and $f$ were defined in eqs. (9) and (8), and

$\tilde{A}:=\frac{\partial L}{\partial x_{i}}y_{i}+\frac{\partial
L}{\partial\dot{x}_{i}}\dot{y}_{i}.$ (24)

Let $\frac{d(...)}{d\epsilon}\Big{|}_{\epsilon\rightarrow 0}:=(...)^{\prime}$.
Then eq. (23) becomes

$\displaystyle\frac{d}{d\epsilon}\left[\frac{df}{dt}\right]\bigg{|}_{\epsilon\rightarrow
0}$ $\displaystyle=\frac{df^{\prime}}{dt}$ (25)
$\displaystyle=\tilde{A}^{\prime}+B^{\prime}f+Bf^{\prime}.$

Eq. (25) is clearly similar to eq. (10), and it can be solved in a similar
fashion. Define

$\frac{dG}{dt}:=\left(\tilde{A}^{\prime}+B^{\prime}f\right)e^{-\int Bd\tau}$
(26)

so that

$\displaystyle f^{\prime}$ $\displaystyle=Ge^{\int Bd\tau}$ (27)
$\displaystyle=\left[\int\left(\tilde{A}^{\prime}+B^{\prime}f\right)e^{-\int
Bd\tau_{1}}d\tau_{2}\right]e^{\int Bd\tau_{3}}$
$\displaystyle=\Bigg{[}\int\Bigg{(}\frac{\partial^{2}L}{\partial x_{i}\partial
x_{j}}y_{i}y_{j}+2\frac{\partial^{2}L}{\partial
x_{i}\partial\dot{x}_{j}}y_{i}\dot{y}_{j}+\frac{\partial^{2}L}{\partial\dot{x}_{i}\dot{x}_{j}}\dot{y}_{i}\dot{y}_{j}+2\frac{\partial^{2}L}{\partial
S\partial\dot{x}_{i}}f\dot{y}_{i}$
$\displaystyle+2\frac{\partial^{2}L}{\partial x_{i}\partial
S}fy_{i}+\frac{\partial^{2}L}{\partial S^{2}}f^{2}\Bigg{)}e^{-\int
Bd\tau_{1}}d\tau_{2}\Bigg{]}e^{\int Bd\tau_{3}}.$

It is tempting, at this point, to think that one can simply set $f=0$, because
the first variation of the action vanishes along a stationary path. We must
remember, however, that this condition only applies at the endpoints of the
system’s trajectory (that is, $f(t_{f})=f(0)=0$, but for some abitrary time
$t$, $f(t)\neq 0$). We therefore need to be careful in handling the terms
proportional to $f$ in eq. (27). At some arbitrary time $t$, eq. (27) can be
written in the form

$f(t)e^{-\int_{0}^{t}B(\tau_{1})d\tau_{1}}=\int_{0}^{t}\left(\frac{\partial
L}{\partial x_{i}}y_{i}+\frac{\partial
L}{\partial\dot{x}_{i}}\dot{y}_{i}\right)e^{-\int_{0}^{\tau_{2}}B(\tau_{1})d\tau_{1}}d\tau_{2}.$
(28)

Integration of the RHS by parts produces

$\displaystyle
f(t)e^{-\int_{0}^{t}B(\tau_{1})d\tau_{1}}=\int_{0}^{t}\bigg{(}\frac{\partial
L}{\partial x_{i}}$ $\displaystyle-\frac{d}{dt}\frac{\partial
L}{\partial\dot{x}_{i}}+\frac{\partial L}{\partial S}\frac{\partial
L}{\partial\dot{x}_{i}}\bigg{)}y_{i}e^{-\int_{0}^{\tau_{2}}B(\tau_{1})d\tau_{1}}d\tau_{2}$
(29)
$\displaystyle+\left[e^{-\int_{0}^{\tau_{2}}B(\tau_{1})d\tau_{1}}\frac{\partial
L}{\partial\dot{x}}y\right]^{t}_{0}.$

Because the Euler-Lagrange equation vanishes at any time $t$, the first term
on the RHS of eq. (29) vanishes. Eq. (29) then reduces to

$f(t)e^{-\int_{0}^{t}B(\tau_{1})d\tau_{1}}=e^{-\int_{0}^{t}B(\tau_{1})d\tau_{1}}\left[\frac{\partial
L}{\partial\dot{x}_{i}}y_{i}\right]\bigg{|}_{\tau_{2}=t}-e^{-\int_{0}^{0}B(\tau_{1})d\tau_{1}}\left[\frac{\partial
L}{\partial\dot{x}_{i}}y_{i}\right]\bigg{|}_{\tau_{2}=0}.$ (30)

Because $y_{i}(0)=0$, it follows that

$f(t)=\frac{\partial L}{\partial\dot{x}_{i}}y_{i},$ (31)

where $\frac{\partial L}{\partial\dot{x}_{i}}y_{i}$ is evaluated at
$\tau_{1}=t$. Insertion of eq. (31) into eq. (27) produces

$\displaystyle f^{\prime}=e^{\int_{0}^{t}Bd\tau_{3}}\int_{0}^{t}$
$\displaystyle\bigg{(}\frac{\partial^{2}L}{\partial\dot{x}_{i}\dot{x}_{j}}\dot{y}_{i}\dot{y}_{j}+2\left[\frac{\partial^{2}L}{\partial
x_{i}\partial\dot{x}_{j}}+\frac{\partial^{2}L}{\partial
S\partial\dot{x}_{i}}\frac{\partial
L}{\partial\dot{x}_{j}}\right]y_{i}\dot{y}_{j}+\frac{\partial^{2}L}{\partial
x_{i}\partial x_{j}}y_{i}y_{j}$ (32)
$\displaystyle+2\frac{\partial^{2}L}{\partial x_{i}\partial S}\frac{\partial
L}{\partial\dot{x}_{j}}y_{j}y_{i}+\frac{\partial^{2}L}{\partial
S^{2}}\frac{\partial L}{\partial\dot{x}_{i}}\frac{\partial
L}{\partial\dot{x}_{j}}y_{i}y_{j}\bigg{)}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}.$

Eq. (32) is the general form of the second variational derivative of the
Herglotz action. Note that, if the Lagrangian does not depend on the action,
then eq. (32) reduces to

$f^{\prime}=\int_{0}^{t}\left(\frac{\partial^{2}L}{\partial\dot{x}_{i}\dot{x}_{j}}\dot{y}_{i}\dot{y}_{j}+2\frac{\partial^{2}L}{\partial
x_{i}\partial\dot{x}_{j}}y_{i}\dot{y}_{j}+\frac{\partial^{2}L}{\partial
x_{i}\partial x_{j}}y_{i}y_{j}\right)d\tau_{1},$ (33)

which is the second variation of the action for a conservative system.

In the next section I show, following Ref. 30, that the second variation of
the Herglotz action vanishes for a particular variation when the system’s
trajectory terminates at a kinetic focus. I then show that the Herglotz action
is a saddle point when the system’s trajectory terminates after a kinetic
focus.

##  III The Herglotz action and kinetic foci

In this section I examine the necessary and sufficient conditions for the
second variation of the Herglotz action to vanish, for the first time, at the
first kinetic focus on a given trajectory. The sufficient condition in this
case is just the same as the sufficient condition for conservative systems
proved in Ref. 30. In that paper the authors showed that, given a stationary
path (call it $x_{0}$), there is a point in time ($t_{\rm KF}$) at which the
second variation of the action along $x_{0}$ vanishes, this being the point at
which a nearby stationary path (call it $x_{1}$) coalesces with $x_{0}$ in the
limit that their initial velocities are made equal. The argument in Ref. 30
makes no reference to the specific form of the action, only requiring that
$\delta^{2}S_{0}>0$ for $t<t_{\rm KF}$. As long as $\delta^{2}S_{0}>0$ for
$t<t_{\rm KF}$, where in this case $S$ is the Herglotz action, then by the
same logic the sufficient condition for $\delta^{2}S$ to vanish is that the
given stationary path ends at $t_{\rm KF}$.

The necessary condition will require a little more work, but this can be done
along the same lines as Ref. 30, from which I adapt the following argument.
Let the action of $x_{0}$ be $S_{0}$, and let the action along $x_{1}$ be
$S_{1}$. In this case $S_{0}$ can be expanded in the series

$S_{0}=S_{1}+\delta S_{1}+\delta^{2}S_{1}+\delta^{3}S_{1}+...$ (34)

and $S_{1}$ can be expanded in the series

$S_{1}=S_{0}+\delta S_{0}+\delta^{2}S_{0}+\delta^{3}S_{0}+...$ (35)

In eqs. (34) and (35) I have defined

$\delta^{(n)}S_{m}=\left(-1\right)^{nm}\frac{\epsilon^{n}}{n!}\left[\frac{d^{n}}{d\epsilon^{n}}S\right]\bigg{|}_{\epsilon\rightarrow
0}$ (36)

with $m=0,1$. The factor of $(-1)^{nm}$ in the definition is due to the fact
that $x_{1,i}=x_{0,i}+\epsilon y_{i}$, while the reference path is
$x_{0,i}=x_{1,i}-\epsilon y_{i}$. That is, to get from $x_{0,i}$ to $x_{1,i}$,
one must vary $x_{0,i}$ in the “positive direction" away from $x_{0,i}$, while
to get from $x_{1,i}$ to $x_{0,i}$ requires a variation in the “negative
direction" away from $x_{1,i}$. This has the effect of inducing a negative
sign in $\delta^{3}S_{1}$, as we will see below (see also Ref. 30 for a
discussion of this point).

If we subtract eq. (34) from eq. (35), we get

$S_{1}-S_{0}=\left(S_{0}-S_{1}\right)+\left(\delta S_{0}-\delta
S_{1}\right)+\left(\delta^{2}S_{0}-\delta^{2}S_{1}\right)+\left(\delta^{3}S_{0}-\delta^{3}S_{1}\right)+...$
(37)

Because both paths are stationary, $\delta S_{0}=\delta S_{1}=0$. Then eq.
(37) can be written in the form

$S_{1}-S_{0}=\frac{1}{2}\left(\delta^{2}S_{0}-\delta^{2}S_{1}\right)+\frac{1}{2}\left(\delta^{3}S_{0}-\delta^{3}S_{1}\right)+...$
(38)

From eq. (32),

$\frac{d^{2}S}{d\epsilon^{2}}\bigg{|}_{\epsilon\rightarrow
0}=e^{\int_{0}^{t}Bd\tau_{3}}\int^{t}_{0}\left(P_{ij}\dot{y}_{i}\dot{y}_{j}+2Q_{ij}y_{i}\dot{y}_{j}+r_{ij}y_{i}y_{j}\right)e^{-\int^{\tau_{2}}_{0}Bd\tau_{1}}d\tau_{2},$
(39)

where I have defined

$P_{ij}:=\frac{\partial^{2}L}{\partial\dot{x}_{i}\partial\dot{x}_{j}},$ (40)
$Q_{ij}:=\frac{\partial^{2}L}{\partial
x_{i}\partial\dot{x}_{j}}+\frac{\partial^{2}L}{\partial
S\partial\dot{x}_{i}}\frac{\partial L}{\partial\dot{x}_{j}},$ (41)

and

$r_{ij}:=\frac{\partial^{2}L}{\partial x_{i}\partial
x_{j}}+2\frac{\partial^{2}L}{\partial x_{i}\partial S}\frac{\partial
L}{\partial\dot{x}_{j}}+\frac{\partial^{2}L}{\partial S^{2}}\frac{\partial
L}{\partial\dot{x}_{i}}\frac{\partial L}{\partial\dot{x}_{j}}.$ (42)

Each of these functions can be expanded around $x_{0}$, so that they take the
forms

$P_{ij}(x_{1})\approx P_{ij}(x_{0})+\frac{\partial P_{ij}(x_{0})}{\partial
x_{k}}\left(x^{(1)}_{k}-x^{(0)}_{k}\right)=P_{ij}(x_{0})+\epsilon\frac{\partial
P_{ij}(x_{0})}{\partial x_{k}}y_{k},$ (43) $Q_{ij}(x_{1})\approx
Q_{ij}(x_{0})+\frac{\partial Q_{ij}(x_{0})}{\partial
x_{k}}\left(x^{(1)}_{k}-x^{(0)}_{k}\right)=Q_{ij}(x_{0})+\epsilon\frac{\partial
Q_{ij}(x_{0})}{\partial x_{k}}y_{k},$ (44)

and

$r_{ij}(x_{1})\approx r_{ij}(x_{0})+\frac{\partial r_{ij}(x_{0})}{\partial
x_{k}}\left(x^{(1)}_{k}-x^{(0)}_{k}\right)=r_{ij}(x_{0})+\epsilon\frac{\partial
r_{ij}(x_{0})}{\partial x_{k}}y_{k}$ (45)

to first order in $\epsilon$. From these expansions, we can see that the
secondary Lagrangian, defined in Ref. 44 as

$2\Omega:=P_{ij}\dot{y}_{i}\dot{y}_{j}+2Q_{ij}y_{i}\dot{y}_{j}+r_{ij}y_{i}y_{j},$
(46)

can be written

$2\Omega(x_{1})\approx 2\Omega(x_{0})+2\frac{\partial\Omega(x_{0})}{\partial
x_{k}}\left(x^{(1)}_{k}-x^{(0)}_{k}\right)=2\Omega(x_{0})+2\epsilon\frac{\partial\Omega(x_{0})}{\partial
x_{k}}y_{k}.$ (47)

Strictly speaking, it isn’t necessary to introduce the secondary Lagrangian. I
have done it here in order to reduce clutter in the following equations. Like
the secondary Lagrangian, the exponential integrals can be expanded to linear
order, giving

$e^{\int B^{(1)}d\tau}\approx e^{\int
B^{(0)}d\tau}e^{\int\epsilon\frac{\partial B(x_{0})}{\partial
x_{k}}y_{k}d\tau}=e^{\int B^{(0)}d\tau}\left(1+\epsilon\int\frac{\partial
B(x_{0})}{\partial x_{k}}y_{k}d\tau\right).$ (48)

Therefore, after a bit of simplification, the first term on the RHS of eq.
(38) becomes

$\displaystyle\frac{1}{2}\left(\delta^{2}S_{0}-\delta^{2}S_{1}\right)$
$\displaystyle=-\frac{\epsilon^{3}}{4}e^{\int
Bd\tau_{3}}\bigg{\\{}\int_{0}^{t}\left[2\frac{\partial\Omega}{\partial
x_{k}}y_{k}-2\Omega\left(\int_{0}^{\tau_{2}}\frac{\partial B}{\partial
x_{k}}y_{k}d\tau_{1}\right)\right]e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}$
(49) $\displaystyle+\int^{t}_{0}\frac{\partial B}{\partial
x_{k}}y_{k}d\tau_{3}\int_{0}^{t}2\Omega
e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}\bigg{\\}}+O(\epsilon^{4}),$

where $\Omega$, $\frac{\partial\Omega}{\partial x_{k}}$, $B$, and
$\frac{\partial B}{\partial x_{k}}$ are evaluated on $x_{0}$.

The key point to keep in mind here is that eq. (49) is $O(\epsilon^{3})$ to
lowest order, which ultimately means (see below) that the second variation of
the action will vanish when the system terminates at a kinetic focus. It also
means that, in calculating $S_{0}-S_{1}$, we can’t ignore the
$\delta^{3}S_{0}-\delta^{3}S_{1}$ term in eq. (38). Calculating the third
functional derivative of $S$ (or, equivalently, the second functional
derivative of $f$) leads to

$\frac{df^{\prime\prime}}{dt}=\tilde{A}^{\prime\prime}+B^{\prime\prime}f+2B^{\prime}f^{\prime}+Bf^{\prime\prime},$
(50)

which is solved by

$f^{\prime\prime}=e^{\int
Bd\tau_{3}}\left[\int_{0}^{t}\left(\tilde{A}^{\prime\prime}+B^{\prime\prime}f+2B^{\prime}f^{\prime}\right)e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}\right].$
(51)

From eq. (31), we have

$f^{\prime\prime}=e^{\int
Bd\tau_{3}}\left[\int_{0}^{t}\left(\tilde{A}^{\prime\prime}+B^{\prime\prime}\frac{\partial
L}{\partial\dot{x}_{i}}y_{i}+2B^{\prime}f^{\prime}\right)e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}\right].$
(52)

As before, one can Taylor expand the RHS of this equation around $x_{0}$. This
is a messier calculation than what I have shown above, so the details are
sketched out in Appendix B. The main point is that, to lowest order, and with
$B$, $\Omega$, $\frac{\partial\Omega}{\partial x_{k}}$, and $\frac{\partial
B}{\partial x_{k}}$ evaluated on $x_{0}$,

$\displaystyle S_{1}-S_{0}$
$\displaystyle=\frac{1}{2}\left(\delta^{2}S_{0}-\delta^{2}S_{1}\right)+\frac{1}{2}\left(\delta^{3}S_{0}-\delta^{3}S_{1}\right)$
(53)
$\displaystyle=\frac{\epsilon^{2}}{4}\left(f^{\prime}_{(0)}-f^{\prime}_{(1)}\right)+\frac{\epsilon^{3}}{12}\left(f^{\prime\prime}_{(0)}+f^{\prime\prime}_{(1)}\right)$
$\displaystyle=-\frac{\epsilon^{3}}{4}\bigg{\\{}e^{\int
Bd\tau_{3}}\int_{0}^{t}\left[2\frac{\partial\Omega}{\partial
x_{k}}y_{k}-2\Omega\int_{0}^{\tau_{2}}\frac{\partial B}{\partial
x_{k}}y_{k}d\tau_{1}\right]e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}$
$\displaystyle+e^{\int Bd\tau_{3}}\int_{0}^{t}\frac{\partial B}{\partial
x_{k}}y_{k}d\tau_{3}\int_{0}^{t}2\Omega
e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}\bigg{\\}}$
$\displaystyle+\frac{\epsilon^{3}}{6}\bigg{\\{}e^{\int
Bd\tau_{3}}\bigg{[}\int_{0}^{t}\bigg{(}\tilde{A}^{\prime\prime}+B^{\prime\prime}\frac{\partial
L}{\partial\dot{x}_{i}}y_{i}+2B^{\prime}f^{\prime}\bigg{)}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}\bigg{]}\bigg{\\}}.$

This equation tells us that $S_{1}-S_{0}\sim O(\epsilon^{3})$ to lowest order
when the stationary path terminates near a kinetic focus (conjugate point).
This has two consequences. First, and as mentioned above, $\delta^{2}S_{0}=0$
when $t_{f}=t_{\rm KF}$. Second, the action of the stationary path $x_{0}$ is
a saddle point when the stationary path terminates a short time after $t_{\rm
KF}$.

That $\delta^{2}S_{0}=0$ when $t_{f}=t_{\rm KF}$ follows from a comparison of
eqs. (38) and (35). Subtraction of $S_{0}$ from both sides of eq. (35)
produces

$S_{1}-S_{0}=\delta^{2}S_{0}+\delta^{3}S_{0}+...$ (54)

Because $S_{1}-S_{0}$ in eq. (38) has the form
$-\frac{\epsilon^{3}}{4}\\{...\\}+\frac{\epsilon^{3}}{6}\\{...\\}$, and
$\delta^{3}S_{0}=-\frac{\epsilon^{3}}{12}\frac{d^{3}S_{0}}{d\epsilon^{3}}$, it
follows that $\delta^{2}S_{0}$ must also be $O(\epsilon^{3})$ if $t_{f}$ is
near $t_{\rm KF}$. The difference $S_{1}-S_{0}$, which is $O(\epsilon^{3})$
when $t_{f}$ is near $t_{\rm KF}$, vanishes when $\epsilon\rightarrow 0$ (that
is, when $x_{1}$ and $x_{0}$ coalesce). Therefore $\delta^{2}S_{0}$, which is
also $O(\epsilon^{3})$, must also vanish when $\epsilon\rightarrow 0$. If the
signs and relative magnitudes of the quantities in the curly brackets are such
that $S_{1}-S_{0}>0$ for a variation that terminates before $t_{\rm KF}$, then
$S_{0}$ will be minimized for $t<t_{\rm KF}$. We have already seen that
$S_{1}-S_{0}=0$ when $\epsilon\rightarrow 0$, which occurs when the varied
path terminates at a kinetic focus. When the varied path terminates after a
kinetic focus, then $\epsilon y_{k}$ will be negative, ensuring that the
varied path $x_{1}$ has an action which is smaller than the action of $x_{0}$
(see Fig. 12 of Ref. 30). $\delta^{2}S_{0}$ will then be negative for a
variation away from $x_{0}$ toward $x_{1}$, but will remain positive for all
other variations. $S_{0}$ is therefore a saddle point when $x_{0}$ terminates
slightly later than $t_{\rm KF}$. As an example, consider the Herglotz
Lagrangian

$L=\frac{1}{2}m\dot{x}^{2}-U(t,x)+g(t)S.$ (55)

Many dissipative systems studied in the literature are governed by Lagrangians
of the form shown in eq. (55),Sloan (2021); Lazo  _et al._ (2018, 2017);
Georgieva and Guenther (2005); Georgieva, Guenther, and Bodurov (2003); Paiva,
Lazo, and Zanchin (2021); Zhang and Tian (2019); Tian and Zhang (2020) as are
the two examples that I study in sections IV.1 and IV.2, so it is worthwhile
to consider this special case. In this case all partial derivatives of $B$
with respect to $x$ vanish, $2\frac{\partial\Omega(x_{0})}{\partial
x}=-U^{\prime\prime\prime}(t,x_{0})y^{3}$, and
$\tilde{A}^{\prime\prime}=-U^{\prime\prime\prime}(t,x_{0})y^{3}$. Then eq.
(38) becomes

$S_{1}-S_{0}=\frac{\epsilon^{3}}{12}e^{\int_{0}^{t}B^{(0)}d\tau_{3}}\int_{0}^{t}U^{\prime\prime\prime}(t,x_{0})y^{3}e^{-\int_{0}^{\tau_{2}}B^{(0)}d\tau_{1}}d\tau_{2}.$
(56)

If $U^{\prime\prime\prime}(t,x)>0$ up to time $t$, then the integrand in eq.
(56) will be positive, and the sign of $S_{1}-S_{0}$ will only change when
$\epsilon y$ changes sign, just beyond $t_{\rm KF}$. Writing $S_{1}-S_{0}$ in
this form not only makes the behavior of the action clearer, but it also makes
it easier to see the relation between the action of a non-conservative, one
dimensional system described by eq. (55) and the action of a conservative one-
dimensional system studied in Ref. 30. $B$ vanishes when $g(t)=0$, in which
case the results of this section reduce to those of Ref. 30 (assuming a time-
independent potential energy function).

To summarize the results of this section, the second variation of the Herglotz
action of a non-conservative system vanishes for the first time when the
trajectory of the system terminates at that trajectory’s first kinetic focus.
The Herglotz action is a saddle point ($\delta^{2}S<0$ for at least one
variation, $\delta^{2}S>0$ for all others) if the trajectory terminates after
the first kinetic focus. These statements generalize the results obtained in
Refs. 30 and 35 for conservative systems, showing that the analysis of
$\delta^{2}S$ for a non-conservative system described by the Herglotz action
can proceed along the same lines as the analysis of $\delta^{2}S$ for a
conservative system, only with somewhat more complicated equations. One
convenient method for calculating the sign of $\delta^{2}S$ for conservative
systems is described in Refs. 36 and 45. In this method, the equation for
$\delta^{2}S$ (eq. 33 here) is transformed into a linear eigenvalue equation,
with the signs of the eigenvalues reflecting the sign of $\delta^{2}S$ at a
particular value of $t_{f}$. I will apply this method to eq. (27) in the
following section, and demonstrate its use with two examples.

##  IV Eigenvalue analysis of second variation of Herglotz action

First, because we’re interested in the sign of $f^{\prime}$, we needn’t worry
about the $e^{\int Bd\tau_{2}}$ term outside the integral. Because the
exponential is positive-definite along the system’s entire trajectory, it will
affect the amplitude of $f^{\prime}$ but not the location of its roots, so it
won’t have any effect on when the sign of $f^{\prime}$ changes (assuming, as I
do throughout this paper, that $B(t)$ is a real-valued function). We can
therefore define a new quantity

$\delta^{2}\mathcal{A}:=f^{\prime}e^{-\int Bd\tau_{3}}$ (57)

and determine how the sign of this quantity changes with $t$. We want to
determine the sign of

$\displaystyle\delta^{2}\mathcal{A}=\int_{0}^{t}\bigg{(}\frac{\partial^{2}L}{\partial\dot{x}_{i}\dot{x}_{j}}\dot{y}_{i}\dot{y}_{j}$
$\displaystyle+2\left[\frac{\partial^{2}L}{\partial
x_{i}\partial\dot{x}_{j}}y_{i}+\frac{\partial^{2}L}{\partial
S\partial\dot{x}_{i}}\frac{\partial
L}{\partial\dot{x}_{j}}y_{i}\right]\dot{y}_{j}+\frac{\partial^{2}L}{\partial
x_{i}\partial x_{j}}y_{i}y_{j}$ (58)
$\displaystyle+2\frac{\partial^{2}L}{\partial x_{i}\partial S}\frac{\partial
L}{\partial\dot{x}_{j}}y_{j}y_{i}+\frac{\partial^{2}L}{\partial
S^{2}}\frac{\partial L}{\partial\dot{x}_{i}}\frac{\partial
L}{\partial\dot{x}_{j}}y_{i}y_{j}\bigg{)}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}.$

Following Refs. 36 and 45, I will define

$R_{ij}:=r_{ij}+BQ_{ji}+\frac{B^{2}}{4}P_{ij}-\frac{1}{2}\frac{dB}{dt}P_{ij}-\frac{B}{2}\frac{dP_{ij}}{dt},$
(59)

and

$\Lambda_{ij}:=-\frac{d}{dt}\left(P_{ij}\frac{d}{dt}+Q_{ji}\right)+\left(Q_{ij}\frac{d}{dt}+R_{ij}\right),$
(60)

where $P_{ij}$, $Q_{ij}$ and $r_{ij}$ are defined in eqs. (40-42). One can
then show (see Appendix A), that eq. (58) takes the form

$\delta^{2}\mathcal{A}=\int_{0}^{t}\eta_{i}\Lambda_{ij}\eta_{j}\hskip
1.42262ptd\tau_{2}$ (61)

where the $\\{\eta_{i}\\}$ functions are defined to be

$\eta_{i}:=y_{i}e^{-\frac{1}{2}\int_{0}^{t}Bd\tau_{1}}.$ (62)

By defining the inner product

$\left(\eta_{i},\eta_{j}\right):=\int\eta_{i}\eta_{j}\hskip 1.42262ptd\tau,$
(63)

and then expanding the $\\{\eta_{i}\\}$ functions in an orthonormal basis
$\\{b_{in}\\}$

$\eta_{i}=\sum_{n}\beta_{n}b_{in},$ (64)

with

$\int b_{n}b_{m}d\tau=\delta_{nm},$ (65)

then eq. (58) takes the form

$\delta^{2}\mathcal{A}=\sum_{n,m}A_{nm}\beta_{n}\beta_{m},$ (66)

where

$A_{nm}=\left(b_{n},\Lambda b_{m}\right)=\int b_{in}\Lambda_{ij}b_{jm}d\tau.$
(67)

If the basis diagonalizes $A_{nm}$, then

$\Lambda_{ij}b_{jn}=\lambda_{n}b_{jn}.$ (68)

We have now reduced the problem of determining the sign of
$\delta^{2}\mathcal{A}$ to a linear eigenvalue problem. If all the eigenvalues
are positive at time $t$, then the action is minimized up to that time. If one
of the eigenvalues ($\lambda_{n}$) vanishes at time $t$, then the action along
the trajectory of the $n^{\rm th}$ variation vanishes, and $t$ is a kinetic
focus. If one of the eigenvalues ($\lambda_{n}$) is negative at time $t$, then
the action is maximized along the trajectory of the $n^{\rm th}$
variation.Papastavridis (1983) Determining the sign of the Herglotz action is
then just a matter of working out the eigenvalue spectrum $\\{\lambda_{n}\\}$
from eq. (68). To illustrate the use of this technique, I will make the
simplifying assumption that the Herglotz Lagrangian takes the form shown in
eq. (55). Then eq. (68) becomes

$\left(-P\frac{d^{2}}{dt^{2}}+R\right)b_{n}=\lambda_{n}b_{n},$ (69)

where

$P=\frac{\partial^{2}L}{\partial\dot{x}^{2}}=m,$ (70)

and

$R=\frac{\partial^{2}L}{\partial
x^{2}}+\frac{1}{4}\left(g^{2}-2\frac{dg}{dt}\right)\frac{\partial^{2}L}{\partial\dot{x}^{2}}=-\frac{\partial^{2}U}{\partial
x^{2}}+\frac{m}{4}\left(g^{2}-2\frac{dg}{dt}\right).$ (71)

Eq. (69) is, as described in Ref. 36, directly analogous to the eigenvalue
equation for a particle in a box with an effective mass
$m=\frac{\hbar^{2}}{2P}$ moving in an effective potential

$\mathcal{U}=\begin{cases}R&\text{if}\ 0<t<t_{f},\\\
\vspace{1mm}\infty&\text{if}\ t=0\hskip 2.84526pt{\rm or}\hskip
2.84526ptt=t_{f},\\\ \end{cases}$ (72)

with $t$ playing the role that $x$ ordinarily plays in time-independent
perturbation theory. Inserting eqs. (70) and (71) into eq. (69), we find

$\left(-\frac{d^{2}}{dt^{2}}-\frac{\partial^{2}\tilde{U}}{\partial
x^{2}}+\frac{g^{2}}{4}-\frac{1}{2}\frac{dg}{dt}\right)b_{n}=\frac{\lambda_{n}}{m}b_{n},$
(73)

where $\tilde{U}:=\frac{U}{m}$. In general, $\tilde{U}$ is a function of the
solution $x(t)$ of the equation of motion given by eq. (17). If this solution
is known, then one can write $\tilde{U}$ explicitly as a function of time, and
then attempt to solve eq. (73). In the special case that
$\frac{\partial^{2}\tilde{U}}{\partial x^{2}}=F(t)+C$, where $F(t)$ is a
function only of time and $C$ is a constant, one can solve the eigenvalue
spectrum of the non-conservative system without solving the equation of motion
beforehand. Eq. (73) can be solved in a number of ways that are familiar from
non-relativistic quantum mechanics, such as perturbation theory, the
variational method, the WKB method, and others.Hussein  _et al._ (1980);
Schumacher, Benjamin and Westmoreland, Michael D. (2010); Shankar, Ramamurti
(1994) In the following sections I will consider two examples: a harmonic
oscillator with time-independent damping, and a harmonic oscillator with time-
dependent damping (and time-dependent frequency).

###  IV.1 Harmonic oscillator with time-independent damping

For an oscillator with time-independent damping, let $g(t)=-a$ and
$\tilde{U}=\frac{1}{2}\omega_{0}^{2}x^{2}$, where $a$ and $\omega_{0}$ are
real-valued constants. In this case eq. (73) reduces to

$\left(-\frac{d^{2}}{dt^{2}}-\omega_{d}^{2}\right)b_{n}=\frac{\lambda_{n}}{m}b_{n},$
(74)

where I have defined $\omega_{d}^{2}:=\omega_{0}^{2}-\frac{a^{2}}{4}$. From
Refs. 30 and 36, we know that the action eigenvalue spectrum of an undamped
harmonic oscillator is

$\frac{\lambda_{n}^{c}}{m}=\left(\frac{n\omega}{2}\right)^{2}-\omega_{0}^{2},$
(75)

where $\omega:=\frac{2\pi}{t}$, so the eigenvalue spectrum of the damped
oscillator must be

$\frac{\lambda_{n}}{m}=\left(\frac{n\omega}{2}\right)^{2}-\omega_{0}^{2}+\frac{a^{2}}{4}.$
(76)

This leads to a familiar result. As is well known, the frequency of a harmonic
oscillator subject to linear damping, with an equation of motion given by

$\ddot{x}+a\dot{x}+\omega_{0}^{2}x=0,$ (77)

is $\omega_{\rm d}=\omega_{0}\sqrt{1-\left(\frac{a}{2\omega_{0}}\right)^{2}}$
(see Ref. 2, for example). If $n=1$, then the second variation of the action
vanishes when

$0=\left(\frac{\omega}{2}\right)^{2}-\omega_{0}^{2}+\frac{a^{2}}{4},$ (78)

or

$\omega=2\omega_{0}\sqrt{1-\left(\frac{a}{2\omega_{0}}\right)^{2}}=2\omega_{d}.$
(79)

From eq. (18), one can show that the first kinetic focus of a damped harmonic
oscillator occurs when

$t_{\rm
KF,1}=\frac{2\pi}{\omega}=\frac{\pi}{\omega_{d}}=\frac{\pi}{\omega_{0}}\frac{1}{\sqrt{1-\left(\frac{a}{2\omega_{0}}\right)^{2}}}=\frac{T_{0}}{2}\frac{1}{\sqrt{1-\left(\frac{a}{2\omega_{0}}\right)^{2}}},$
(80)

which is just the half-period of the damped oscillator ($T_{0}$ is the
undamped period).555This result was also obtained by G. Leitmann in Ref. 15,
with an action of the form $S=\int
e^{2ct}\left[\dot{x}^{2}-\frac{1}{2}\omega_{0}^{2}x^{2}\right]dt$, where my
$a$ equals his $2c$. Eq. (80) also tells us that the time to reach the first
kinetic focus diverges as the oscillator approaches critical damping (that is,
when $a\rightarrow 2\omega_{0}$). In this case (or in the case of strong
damping, for which $a>2\omega_{0}$) the oscillations decay on a time scale
that is shorter than $T_{0}/2$, so there are no kinetic foci, and the action
is minimized along the entire trajectory of the oscillator.

###  IV.2 Harmonic oscillator with time-independent damping

Consider an equation of motion of the form

$\ddot{x}+\gamma(t)\dot{x}+\omega^{2}(t)x=0,$ (81)

where $\gamma(t)$ is the dissipation function and the frequency depends on
time. This equation describes a time-dependent harmonic oscillator (TDHO), a
system that has been used to model the motion of a non-relativistic charged
particle in a time-dependent magnetic field Lewis (1967) as well as the normal
mode dynamics of a scalar field in an expanding universeRobles-Pérez (2017);
Finelli, Vacca, and Venturi (1998) (see also Ref. 52).

Eq. (81) can be derived from the Herglotz Lagrangian of an oscillator with
time-dependent frequency and unit mass

$L=\frac{1}{2}\dot{x}^{2}-\frac{1}{2}\omega^{2}(t)x^{2}+g(t)S,$ (82)

where, in this case, $g(t)=-\gamma(t)$. Eq. (73) tells us that the action
eigenvalue equation for the TDHO is

$\left(-\frac{d^{2}}{dt^{2}}-\omega^{2}(t)+\frac{1}{4}\gamma^{2}(t)+\frac{1}{2}\dot{\gamma}(t)\right)b_{n}=\lambda_{n}b_{n}.$
(83)

To simplify the problem, I will assume that the damping function and the
frequency vary slowly over the time interval $0\leq t\leq t_{f}$. I will also
assume that these functions can be Taylor expanded around $t=0$, so that

$\omega(t)=\omega(0)+t\frac{d\omega}{dt}\bigg{|}_{t=0}+O(t^{2}),$ (84)

and

$\gamma(t)=\gamma(0)+t\frac{d\gamma}{dt}\bigg{|}_{t=0}+O(t^{2}).$ (85)

Define $\omega_{0}:=\omega(0)$,
$\omega_{0}^{\prime}=\frac{d\omega}{dt}\big{|}_{t=0}$,
$\gamma_{0}:=\gamma(0)$, and
$\gamma_{0}^{\prime}=\frac{d\gamma}{dt}\big{|}_{t=0}$. Then eq. (83) takes the
form

$\left(-\frac{d^{2}}{dt^{2}}-\omega_{0}^{2}+\frac{\gamma_{0}^{2}}{4}+\frac{\gamma_{0}^{\prime}}{2}-2\left(\omega_{0}\omega_{0}^{\prime}-\frac{\gamma_{0}\gamma_{0}^{\prime}}{4}\right)t-\left(\omega^{\prime
2}_{0}-\frac{\gamma^{\prime
2}_{0}}{4}\right)t^{2}\right)b_{n}=\lambda_{n}b_{n}.$ (86)

Here we have an effective Schrödinger equation for a particle of unit mass
moving in an effective potential

$\mathcal{U}=\begin{cases}-\omega_{0}^{2}+\frac{\gamma_{0}^{2}}{4}+\frac{\gamma_{0}^{\prime}}{2}-2\left(\omega_{0}\omega_{0}^{\prime}-\frac{\gamma_{0}\gamma_{0}^{\prime}}{4}\right)t-\left(\omega^{\prime
2}_{0}-\frac{\gamma^{\prime 2}_{0}}{4}\right)t^{2}&\text{if}\ 0<t<t_{f},\\\
\vspace{1mm}\infty&\text{if}\ t=0\hskip 2.84526pt{\rm or}\hskip
2.84526ptt=t_{f}.\\\ \end{cases}$ (87)

We can see that, depending on the initial values of the frequency and the
damping function (as well as the initial values of their first derivatives),
the depth of the effective potential can increase with time. Therefore, even
if the system is critically damped at first (so that
$\omega_{0}=2\gamma_{0}$), it may eventually develop a kinetic focus. This
happens because, for a sufficiently deep effective potential, one or more of
the action eigenvalues can be $\leq 0$ (the analogous situation in quantum
mechanics would be to a particle in a potential well that is deep enough to
admit energy eigenvalues that are $\leq 0$).

Figure 2: Lowest action eigenvalue, normalized to $\omega_{0}^{2}$, versus $\omega_{0}t_{f}$, where $t_{f}$ is the time at which the stationary path of the TDHO terminates. The solid pink curve corresponds to a time-independent damped harmonic oscillator, and the other curves correspond to the trial functions shown in Table 3. TDHO trajectories | $v(t=0)$
---|---
$x_{0}(t)$ | 1
$x_{1}(t)$ | 1.05
$x_{2}(t)$ | 1.10
$x_{3}(t)$ | 1.15
$x_{4}(t)$ | 1.20
Table 2: Trajectories used to generate Figs. 3 and 4. The left column lists the trajectories with their corresponding label subscripts, and the right column lists the initial values of their velocities. Variational parameters | Trial function | Legend label
---|---|---
$\alpha$ | $N\left(1+\alpha t\right){\rm sin}\left(\frac{\pi t}{t_{f}}\right)$ | $\psi_{1}$
$\alpha$ | $Ne^{\alpha t}{\rm sin}\left(\frac{\pi t}{t_{f}}\right)$ | $\psi_{2}$
$\alpha$ | $Ne^{\alpha t^{2}}{\rm sin}\left(\frac{\pi t}{t_{f}}\right)$ | $\psi_{3}$
$\alpha$, $\beta$ | $N\left(\alpha+\beta t\right){\rm sin}\left(\frac{\pi t}{t_{f}}\right)$ | $\psi_{4}$
$\alpha$, $\beta$, $\gamma$ | $N\left(\alpha+\beta t+\gamma t^{2}\right){\rm sin}\left(\frac{\pi t}{t_{f}}\right)$ | $\psi_{5}$
$\alpha$, $\beta$, $\gamma$, $\delta$ | $N\left(\alpha+\beta t+\gamma t^{2}+\delta t^{3}\right){\rm sin}\left(\frac{\pi t}{t_{f}}\right)$ | $\psi_{6}$
Table 3: Trial functions for the variational calculation of the lowest action
eigenvalue of the TDHO. $N$ is a normalization constant and $t_{f}$ is the
time at which the stationary path terminates.

To investigate this possibility, I will estimate the lowest eigenvalue of the
action spectrum, as a function of the final time $t_{f}$, via the variational
method. In standard quantum mechanics, it is possible to estimate the ground
state energy $E_{0}$ of a system with a given Hamiltonian $H$ by calculating
the quantity

$E(\alpha)=\frac{\bra{\psi_{\alpha}}H\ket{\psi_{\alpha}}}{\bra{\psi_{\alpha}}\ket{\psi_{\alpha}}},$
(88)

and demanding that $\alpha$ take its minimum value (in this equation
$\ket{\psi_{\alpha}}$ is a state that depends on the variational parameter
$\alpha$). Then $E(\alpha_{\rm min})$ provides an upper bound on the ground
state energy such that $E(\alpha_{\rm min})\geq E_{0}$.Schumacher, Benjamin
and Westmoreland, Michael D. (2010); Shankar, Ramamurti (1994) The accuracy of
the upper bound $E(\alpha_{\rm min})$ depends on how closely the trial
wavefunction $\psi_{\alpha}$ approximates the ground state wavefunction
$\psi_{0}$ (although it needn’t be too close; even a relatively crude trial
wavefunction can still give a reasonable approximation to the ground state
energy Schumacher, Benjamin and Westmoreland, Michael D. (2010); Shankar,
Ramamurti (1994)). In our case, we want

$\lambda_{1}(\alpha_{\rm
min})=\frac{\left(\psi_{\alpha},\mathcal{H}\psi_{\alpha}\right)}{\left(\psi_{\alpha},\psi_{\alpha}\right)}\bigg{|}_{\alpha=\alpha_{\rm
min}},$ (89)

where the effective Hamiltonian is

$\mathcal{H}=-\frac{d^{2}}{dt^{2}}+\mathcal{U}.$ (90)

To compute $\lambda_{1}(\alpha_{\rm min})$, I wrote a Python script to
minimize $\lambda_{1}(\alpha)$ by way of the minimize routine built into
SciPy.666When I called this routine, I used the Sequential Least Squares
Programming method. Documentation:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
I used several trial functions, shown in Table 3. The ${\rm
sin}\left(\frac{\pi t}{t_{f}}\right)$ in each trial function is to ensure that
the function vanishes at $t=0$ and $t=t_{f}$. For the single parameter
functions, I wanted to ensure that the trial function also decays rapidly in
the direction of $t=0$, hence the choice of a linear, exponential, or Gaussian
function multiplying the sine. The multi-parameter functions are intended as a
convergence test. The main difficulty with the variational method in quantum
mechanics is that one can never be completely sure of how well it approximates
the true ground state energy of a given system, and the same is true of the
lowest order action eigenvalue here. Generally speaking, if one finds that the
variational estimate does not change significantly as one uses increasingly
complicated trial functions, then one can be fairly confident that the
estimate is reliable.Schumacher, Benjamin and Westmoreland, Michael D. (2010);
Shankar, Ramamurti (1994) The results of this analysis are shown in Fig. 2. To
plot the results in a scale-independent fashion, I divided both sides of eq.
(86) by $\omega_{0}^{2}$ to obtain

$\displaystyle\Bigg{[}-\frac{d^{2}}{d\left(\omega_{0}t\right)^{2}}-1+\frac{\gamma_{0}^{2}}{4\omega_{0}^{2}}+\frac{\gamma_{0}^{\prime}}{2\omega_{0}^{2}}$
$\displaystyle-2\left(\frac{\omega_{0}^{\prime}}{\omega_{0}^{2}}-\frac{\gamma_{0}\gamma_{0}^{\prime}}{4\omega_{0}^{3}}\right)\left(\omega_{0}t\right)$
(91) $\displaystyle-\left(\frac{\omega^{\prime
2}_{0}}{\omega_{0}^{4}}-\frac{\gamma^{\prime
2}_{0}}{4\omega_{0}^{4}}\right)\left(\omega_{0}t\right)^{2}\Bigg{]}b_{n}=\frac{\lambda_{n}}{\omega_{0}^{2}}b_{n}.$

and used the effective Hamiltonian shown on the LHS of eq. (91) in the
variational calculation of eq. (89). In Fig. 2, I plot
$\frac{\lambda_{1}}{\omega_{0}^{2}}$ versus $\omega_{0}t_{f}$, after (somewhat
arbitrarily) setting $\gamma_{0}=2\omega_{0}$,
$\omega_{0}^{\prime}=\omega_{0}^{2}/2$,
$\gamma_{0}^{\prime}=\omega_{0}^{2}/10$, and making the replacements
$t\rightarrow\omega_{0}t$ and $t_{f}\rightarrow\omega_{0}t_{f}$ in the trial
functions of Table 3. The solid pink curve corresponds to a critically damped,
time-independent harmonic oscillator. Each of the other curves corresponds to
one of the trial functions shown in Table 3, labeled accordingly. From the
figure, we can see that value of $\frac{\lambda_{1}}{\omega_{0}^{2}}$ for a
TDHO differs noticeably from that of the time-independent harmonic oscillator,
even down to $\omega_{0}t_{f}\sim 1$. Near $\omega_{0}t_{f}\approx 2.5$, the
eigenvalue curve passes through zero, signaling the appearance of a kinetic
focus. We can also see that all of the trial functions are in good agreement
up to $\omega_{0}t_{f}\approx 3.5$, showing that the prediction of a kinetic
focus at $\omega_{0}t_{f}\approx 2.5$ is robust, given the set of parameters I
have chosen.

Figure 3: Numerically obtained trajectories of a critically damped TDHO with
slightly different initial velocities. See Table 2 for a list of these initial
velocities. Here the displacement $x(t)$ is reported in arbitrary units.

Figure 4: Trajectories of a critically damped time-independent harmonic
oscillator with slightly different initial velocities. See Table 2 for a list
of these initial velocities. Here the displacement $x(t)$ is reported in
arbitrary units.

To check this method, I numerically solved the equation of motion of the TDHO
(eq. 81) using the differential equation solver odeint packaged with
SciPy,777https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html
and plotted the results in Fig. (3). Each curve has an initial displacement of
$x(t=0)=0$ and an initial velocity given in Table 2, where

$v:=\frac{dx}{d\left(\omega_{0}t\right)}.$ (92)

From Fig. 3, we can see that trajectories having different initial velocities
intersect for the first time around $\omega_{0}t\approx 2.5$, in agreement
with what Fig. 2 predicts.888Like the undamped harmonic oscillator, the TDHO
is a special case in that trajectories having arbitrarily different initial
velocities will intersect at the kinetic focus; it is not necessary for the
difference between the initial velocities of two trajectories to be small. If
one were to solve instead the equation of motion of a critically damped, time-
independent harmonic oscillator, one would find the results shown in Fig. 4.
In this case, trajectories with different initial velocities do not intersect,
again in agreement with the prediction of Fig. 2.

These findings demonstrate that the dynamical stability of a non-conservative
system can, under some circumstances, be qualitatively different from the
dynamical stability of a conservative system, and that the Herglotz action is
amenable to the same kind of stability analysis as the standard conservative
action used in Ref. 36.

##  V Conclusion

The sign of the second variation of the action, for a conservative system, can
be used as an indicator of that system’s dynamical stability. As we have seen
(and as shown in Refs. 30; 32; 35, for example), the second variation changes
sign when a conservative system passes through a kinetic focus, and it is the
existence of these kinetic foci that indicate dynamical stability.

In this paper, I have generalized the analyses of Refs. 30 and 36 by showing
that the second variation of the Herglotz action for a non-conservative system
also changes sign when that system passes through a kinetic focus. This means
that the second variation of the Herglotz action can be used to infer the
dynamical stability properties of the class of non-conservative systems that
can be described by the Herglotz action. Further, the eigenvalue analysis of
the action for conservative systems described in Refs. 36 and 45 can also be
applied to non-conservative systems. By doing this, I was able to locate the
first kinetic focus of a TDHO, and I showed that the dynamical stability of
this system, under a particular set of assumptions, is qualitatively different
from the dynamical stability of a time-independent harmonic oscillator. It
would be interesting to extend these results to more physically motivated
applications of the TDHO, like those mentioned in Sec. IV.2.

Another way to extend these results would be to derive the semiclassical
propagator, in the path integral formalism, for a quantum system that has
dissipation and whose classical action is the Herglotz action. As I mentioned
in the Introduction, the number of kinetic foci along a stationary (classical)
path of a conservative quantum system determines the phase of the
semiclassical propagator for that system. It would therefore be interesting to
use the formalism developed here to investigate how dissipation affects the
phase of the semiclassical propagator for a system whose dynamics can be
derived from the Herglotz action. I leave this for future work.

###### Acknowledgements.

I would like to thank Joel Meyers (SMU) for his detailed and insightful
feedback at various stages of this project. Computational resources for this
project were provided by SMU’s Center for Research Computing, and I used the
following Python packages: Matplotlib,Hunter (2007) NumPy,Harris  _et al._
(2020) and SciPy.Virtanen  _et al._ (2020) I was financially supported, in
part, by DOE grant DE-SC0010129.

##  Appendix A Derivation of the action eigenvalue equation

If we integrate the first term in parentheses on the RHS of eq. (58), we get

$\displaystyle\int_{0}^{t}\dot{y}_{i}P_{ij}\dot{y}_{j}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}$
(93)
$\displaystyle=\int_{0}^{t}y_{i}\left[-\frac{d}{dt}\left(P_{ij}\dot{y}_{j}\right)+BP_{ij}\dot{y}_{j}\right]e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}$
$\displaystyle=\int_{0}^{t}y_{i}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}\left[-\frac{d}{dt}\left(P_{ij}\dot{y}_{j}\right)+BP_{ij}\dot{y}_{j}\right]e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2},$

and if we integrate the second term in parentheses on the RHS of eq. (58), we
get

$\displaystyle\int_{0}^{t}2Q_{ij}y_{i}\dot{y}_{j}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}$
(94)
$\displaystyle=\int_{0}^{t}y_{i}Q_{ij}\dot{y}_{j}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}+\int_{0}^{t}y_{j}Q_{ji}\frac{dy_{i}}{dt}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}$
$\displaystyle=\int_{0}^{t}y_{i}Q_{ij}\dot{y}_{j}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}-\int_{0}^{t}y_{j}\frac{d}{dt}\left(Q_{ji}y_{i}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}\right)d\tau_{2}$
$\displaystyle=\int_{0}^{t}y_{i}Q_{ij}\dot{y}_{j}e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}$
$\displaystyle-\int_{0}^{t}\left(y_{j}\dot{Q}_{ji}y_{i}+y_{j}Q_{ji}\dot{y}_{i}-y_{j}BQ_{ji}y_{i}\right)e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}$
$\displaystyle=\int_{0}^{t}y_{i}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}\left(-\dot{Q}_{ji}+BQ_{ji}\right)y_{j}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}.$

(recall the definitions of $P_{ij}$ and $Q_{ij}$ from eqs. 40 and 41). On the
other hand, if we insert eqs. (59), (60), and (62) directly into eq. (61), we
find

$\displaystyle\delta^{2}\mathcal{A}$
$\displaystyle=\int_{0}^{t}y_{i}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}\left[-\frac{d}{dt}\left(P_{ij}\frac{d}{dt}+Q_{ji}\right)+\left(Q_{ij}\frac{d}{dt}+R_{ij}\right)\right]y_{j}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}$
(95)
$\displaystyle=\int_{0}^{t}y_{i}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}\bigg{\\{}\bigg{[}-\frac{d}{dt}\left(P_{ij}\dot{y}_{j}\right)+BP_{ij}\dot{y}_{j}$
$\displaystyle+\left(\frac{B}{2}\dot{P}_{ij}-\frac{B^{2}}{4}P_{ij}+\frac{1}{2}\frac{dB}{dt}P_{ij}\right)y_{j}\bigg{]}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}-\frac{dQ_{ji}}{dt}y_{j}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}$
$\displaystyle+\bigg{(}\frac{\partial^{2}L}{\partial x_{i}\partial
x_{j}}+2\frac{\partial^{2}L}{\partial x_{i}\partial S}\frac{\partial
L}{\partial\dot{x}_{j}}+\frac{\partial^{2}L}{\partial S^{2}}\frac{\partial
L}{\partial\dot{x}_{i}}\frac{\partial L}{\partial\dot{x}_{j}}$
$\displaystyle+BQ_{ji}+\frac{B^{2}}{4}P_{ij}-\frac{1}{2}\frac{dB}{dt}P_{ij}-\frac{B}{2}\dot{P}_{ij}\bigg{)}y_{j}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}\bigg{\\}}d\tau_{2}.$

After canceling terms, this reduces to

$\displaystyle\delta^{2}\mathcal{A}$
$\displaystyle=\int_{0}^{t}y_{i}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}\bigg{\\{}\bigg{[}-\frac{d}{dt}\left(P_{ij}\dot{y}_{j}\right)+BP_{ij}\dot{y}_{j}\bigg{]}$
(96) $\displaystyle+\bigg{(}\frac{\partial^{2}L}{\partial x_{i}\partial
x_{j}}+2\frac{\partial^{2}L}{\partial x_{i}\partial S}\frac{\partial
L}{\partial\dot{x}_{j}}+\frac{\partial^{2}L}{\partial S^{2}}\frac{\partial
L}{\partial\dot{x}_{i}}\frac{\partial
L}{\partial\dot{x}_{j}}+BQ_{ji}-\frac{dQ_{ji}}{dt}\bigg{)}y_{j}\bigg{\\}}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2},$

which is equal to the sum of eqs. (93), (94), and

$\int_{0}^{t}y_{i}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}\bigg{(}\frac{\partial^{2}L}{\partial
x_{i}\partial x_{j}}+2\frac{\partial^{2}L}{\partial x_{i}\partial
S}\frac{\partial L}{\partial\dot{x}_{j}}+\frac{\partial^{2}L}{\partial
S^{2}}\frac{\partial L}{\partial\dot{x}_{i}}\frac{\partial
L}{\partial\dot{x}_{j}}\bigg{)}y_{j}e^{-\frac{1}{2}\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}.$
(97)

The RHS of eq. (58) is also equal to expression (97) plus eqs. (93) and (94),
so it must equal eq. (96). Therefore $\delta^{2}\mathcal{A}$ can be written in
the form of eq. (61).

##  Appendix B Third variation of the Herglotz action

$\tilde{A}^{\prime\prime}$ consists of various combinations of third
derivative of $L$, of the form

$\tilde{A}^{\prime\prime}\sim\frac{\partial^{3}L}{\partial X_{i}\partial
X_{j}\partial X_{k}}Y_{i}Y_{j}Y_{k},$ (98)

where $X_{i}=x_{i},\dot{x}_{i}$, or $S$, and $Y_{i}=y_{i},\dot{y}_{i}$, or
$f$. To lowest order

$\frac{\partial^{3}L}{\partial X_{i}\partial X_{j}\partial
X_{k}}\left(x_{1}\right)\approx\frac{\partial^{3}L}{\partial X_{i}\partial
X_{j}\partial X_{k}}\left(x_{0}\right)+\epsilon\frac{\partial^{4}L}{\partial
X_{i}\partial X_{j}\partial X_{k}\partial x_{l}}\left(x_{0}\right)y_{l}.$ (99)

In the Taylor expansion of the action, $\delta^{3}S$ is associated with a
factor of $\frac{\epsilon^{3}}{6}$. Therefore, to lowest order, the second
term on the RHS of eq. (99) can be dropped in eq. (52). In a similar fashion,

$B^{\prime\prime}\sim\frac{\partial^{3}L}{\partial X_{i}\partial X_{j}\partial
S}Y_{i}Y_{j},$ (100)

and

$\frac{\partial^{3}L}{\partial X_{i}\partial X_{j}\partial
S}\left(x_{1}\right)\approx\frac{\partial^{3}L}{\partial X_{i}\partial
X_{j}\partial S}\left(x_{0}\right)+\epsilon\frac{\partial^{4}L}{\partial
X_{i}\partial X_{j}\partial S\partial x_{k}}y_{k}.$ (101)

Also,

$\frac{\partial L}{\partial\dot{x}_{i}}\left(x_{1}\right)\approx\frac{\partial
L}{\partial\dot{x}_{i}}\left(x_{0}\right)+\epsilon\frac{\partial
L}{\partial\dot{x}_{i}\partial x_{j}}y_{j}$ (102)

As before, we can drop $O(\epsilon)$ terms from the product
$B^{\prime\prime}\frac{\partial L}{\partial\dot{x}_{i}}$. Finally,

$B^{\prime}\sim\frac{\partial^{2}L}{\partial X_{i}\partial S}Y_{i}$ (103)

and

$\frac{\partial^{2}L}{\partial X_{i}\partial
S}\left(x_{1}\right)\approx\frac{\partial^{2}L}{\partial X_{i}\partial
S}\left(x_{0}\right)+\epsilon\frac{\partial^{3}L}{\partial X_{i}\partial
S\partial x_{j}}y_{j},$ (104)

and we know that

$\displaystyle f^{\prime}$ $\displaystyle=e^{-\int
B^{(0)}d\tau_{3}}\int_{0}^{t}\left[2\Omega(x_{0})\right]e^{-\int_{0}^{\tau_{2}}B^{(0)}d\tau_{1}}d\tau_{2}+O(\epsilon)$
(105) $\displaystyle=:f^{\prime}_{(0)}+O(\epsilon).$

Therefore

$f^{\prime\prime}(x_{1})=e^{\int
Bd\tau_{3}}\left[\int_{0}^{t}\left(\tilde{A}^{\prime\prime}+B^{\prime\prime}\frac{\partial
L}{\partial\dot{x}_{i}}y_{i}+2B^{\prime}f^{\prime}\right)e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}\right]+O(\epsilon)$
(106)

and

$f^{\prime\prime}(x_{0})=e^{\int
Bd\tau_{3}}\left[\int_{0}^{t}\left(\tilde{A}^{\prime\prime}+B^{\prime\prime}\frac{\partial
L}{\partial\dot{x}_{i}}y_{i}+2B^{\prime}f^{\prime}\right)e^{-\int_{0}^{\tau_{2}}Bd\tau_{1}}d\tau_{2}\right],$
(107)

(with $B$, $\tilde{A}^{\prime\prime}$, $B^{\prime}$, $B^{\prime\prime}$,
$f^{\prime}$, and $\frac{\partial L}{\partial\dot{x}_{i}}$ in both equations
evaluated on $x_{0}$).

## References

  * Goldstein, Poole, and Safko (2002) H. Goldstein, C. Poole, and J. Safko,  _Classical Mechanics, $3^{rd}$ ed._ (Addison Wesley, San Francisco, CA, 2002).
  * Fowles, Grant R. and Cassiday, George L. (2005) Fowles, Grant R. and Cassiday, George L.,  _Analytical Mechanics, $7^{th}$ ed._ (Thomson Brooks/Cole, Belmont, CA, 2005).
  * Fitzgerald, Massoudi, and Vaidya (2021) K. Fitzgerald, M. Massoudi, and A. Vaidya, “On the modified least action principle with dissipation,” European Journal of Mechanics B / Fluids 89, 301–311 (2021).
  * Riewe (1996) F. Riewe, “Nonconservative lagrangian and hamiltonian mechanics,” Physical Review E 53, 1890 (1996).
  * Galley (2013) C. R. Galley, “Classical Mechanics of Nonconservative Systems,” Physical Review Letters 110, 174301 (2013), arXiv:1210.2745 [gr-qc] .
  * Lemos (1981) N. A. Lemos, “Note on the Langrangian description of dissipative systems,” American Journal of Physics 49, 1181–1182 (1981).
  * Musielak (2008) Z. E. Musielak, “Standard and non-standard Lagrangians for dissipative dynamical systems with variable coefficients,” Journal of Physics A: Mathematical and General 41, 055205 (2008).
  * Denman (1966) H. H. Denman, “On Linear Friction in Lagrange’s Equation,” American Journal of Physics 34, 1147–1149 (1966).
  * Wang and Wang (2012) Q. A. Wang and R. Wang, “Is it possible to formulate least action principle for dissipative systems?” arXiv e-prints , arXiv:1201.6309 (2012), arXiv:1201.6309, arXiv:1201.6309 [cond-mat.stat-mech] .
  * Lin and Wang (2013) T. Lin and Q. A. Wang, “The Extrema of an Action Principle for Dissipative Mechanical Systems,” Journal of Applied Mechanics 81, 031002 (2013), arXiv:1310.0455 [physics.class-ph] .
  * Udwadia, Leitmann, and Cho (2011) F. E. Udwadia, G. Leitmann, and H. Cho, “Some Further Remarks on Hamilton’s Principle,” Journal of Applied Mechanics 78, 011014 (2011).
  * Kobussen (1979) J. A. Kobussen, “Some Comments on the Lagrangian Formalism for Systems with General Velocity-Dependent Forces,” Acta physica Austriaca 51, 293–309 (1979).
  * Bahar and Kwatny (1981) L. Y. Bahar and H. G. Kwatny, “Generalized Lagrangian and conservation law for the damped harmonic oscillator,” American Journal of Physics 49, 1062–1065 (1981).
  * Langley (2021) R. S. Langley, “A generalization of the Lagrange-Hamilton formalism with application to non-conservative systems and the quantum to classical transition,” Journal of Mathematical Physics 62, 033503 (2021).
  * Leitmann (1963) G. Leitmann, “Some Remarks on Hamilton’s Principle,” Journal of Applied Mechanics 30, 623–625 (1963).
  * Gray, Karl, and Novikov (2004) C. G. Gray, G. Karl, and V. A. Novikov, “Progress in classical and quantum variational principles,” Reports on Progress in Physics 67, 159–208 (2004), arXiv:physics/0312071 [physics.class-ph] .
  * Note (1) See, for example, Refs. 3; 10; 59; 7; 4; 60; 61; 6; 62; 63; 64 and references therein.
  * Lazo  _et al._ (2018) M. J. Lazo, J. Paiva, J. T. S. Amaral, and G. S. F. Frederico, “An action principle for action-dependent Lagrangians: Toward an action principle to non-conservative systems,” Journal of Mathematical Physics 59, 032902 (2018), arXiv:1803.08308 [math-ph] .
  * Guenther, R. B., Guenther, C. M., and Gottsch, J. A. (1996) Guenther, R. B., Guenther, C. M., and Gottsch, J. A.,  _The Herglotz Lectures on Contact Transformations and Hamiltonian Systems, Volume 1 of Lecture Notes in Nonlinear Analysis_ (Juliusz Schauder Center for Nonlinear Studies, Nicholas Copernicus University, Toruń, Poland, 1996).
  * Georgieva, Guenther, and Bodurov (2003) B. Georgieva, R. Guenther, and T. Bodurov, “Generalized variational principle of Herglotz for several independent variables. First Noether-type theorem,” Journal of Mathematical Physics 44, 3911–3927 (2003).
  * Lazo  _et al._ (2017) M. J. Lazo, J. Paiva, J. T. S. Amaral, and G. S. F. Frederico, “Action principle for action-dependent lagrangians toward nonconservative gravity: Accelerating universe without dark energy,” Physical Review D 95, 101501 (2017).
  * Paiva, Lazo, and Zanchin (2021) J. A. P. Paiva, M. J. Lazo, and V. T. Zanchin, “Generalized nonconservative gravitational field equations from Herglotz action principle,” arXiv e-prints , arXiv:2108.02902 (2021), arXiv:2108.02902, arXiv:2108.02902 [gr-qc] .
  * Georgieva and Guenther (2005) B. Georgieva and R. Guenther, “Second Noether-type theorem for the generalized variational principle of Herglotz,” Topological Methods in Nonlinear Analysis, Journal of the Juliusz Schauder Center 26, 307–314 (2005).
  * Zhang and Tian (2019) Y. Zhang and X. Tian, “Conservation laws of nonconservative nonholonomic system based on Herglotz variational problem,” Physics Letters A 383, 691–696 (2019).
  * Tian and Zhang (2020) X. Tian and Y. Zhang, “Adiabatic Invariants of Herglotz Type for Perturbed Nonconservative Lagrangian Systems,” Theoretical and Mathematical Physics 202, 126–135 (2020).
  * Sloan (2021) D. Sloan, “New action for cosmology,” Physical Review D 103, 043524 (2021), arXiv:2010.07329 [gr-qc] .
  * Donchev (2014) V. Donchev, “Variational symmetries, conserved quantities and identities for several equations of mathematical physics,” Journal of Mathematical Physics 55, 032901 (2014).
  * Gaset  _et al._ (2020) J. Gaset, X. Gràcia, M. C. Muñoz-Lecanda, X. Rivas, and N. Román-Roy, “New contributions to the Hamiltonian and Lagrangian contact formalisms for dissipative mechanical systems and their symmetries,” International Journal of Geometric Methods in Modern Physics 17, 2050090-39 (2020), arXiv:1907.02947 [math-ph] .
  * Caramês  _et al._ (2018) T. R. P. Caramês, H. Velten, J. C. Fabris, and M. J. Lazo, “Dark energy with zero pressure: Accelerated expansion and large scale structure in action-dependent Lagrangian theories,” Physical Review D 98, 103501 (2018), arXiv:1808.02798 [gr-qc] .
  * Gray and Taylor (2007) C. G. Gray and E. F. Taylor, “When action is not least,” American Journal of Physics 75, 434–458 (2007).
  * Papastavridis (1982) J. G. Papastavridis, “A direct variational method for non-conservative system stability,” Journal of Sound and Vibration 80, 447–459 (1982).
  * Papastavridis (1983) J. G. Papastavridis, “An eigenvalue criterion for the study of the Hamilton action’s extremality,” Mechanics Research Communications 10, 171–179 (1983).
  * Papastavridis (1985) J. G. Papastavridis, “An Action-Based Method for the Kinetic Stability of Potential But Nonconservative Systems,” Journal of Applied Mechanics 52, 731 (1985).
  * Note (2) See Ref. 45. Technically the phase is equal to $-\nu\pi/2$, where $\nu$ is the number of negative eigenvalues of the operator $\Lambda_{ij}$, which I define in eq. (60).
  * Gray and Poisson (2011) C. G. Gray and E. Poisson, “When action is not least for orbits in general relativity,” American Journal of Physics 79, 43–56 (2011).
  * Hussein  _et al._ (1980) M. S. Hussein, J. G. Pereira, V. Stojanoff, and H. Takai, “The sufficient condition for an extremum in the classical action integral as an eigenvalue problem,” American Journal of Physics 48, 767–770 (1980).
  * Papastavridis and Chen (1986) J. G. Papastavridis and G. Chen, “The principle of “least” action in non-linear and/or non-conservative oscillations,” Journal of Sound and Vibration 109, 225–235 (1986).
  * Note (3) An explicit calculation of the second variation of the Herglotz action, similar to what I derive in Sec. II.2, can be found in Ref. 65. As far as I can tell, however, the authors of that paper appear to have been concerned only with proving that the solution of the Herglotz variational problem minimizes the Herglotz action, and not with investigating the sign change of the second variation. Additionally, these authors do not appear to make the assumption, as I do, that the second variation of the Herglotz action vanishes at the initial time $t_{i}$.
  * Vermeeren, Bravetti, and Seri (2019) M. Vermeeren, A. Bravetti, and M. Seri, “Contact variational integrators,” Journal of Physics A: Mathematical and Theoretical 52, 445206 (2019), arXiv:1902.00436 [math.NA] .
  * Gaset and Marín-Salvador (2021) J. Gaset and A. Marín-Salvador, “Application of Herglotz’s Variational Principle to Electromagnetic Systems with Dissipation,” arXiv e-prints , arXiv:2108.07542 (2021), arXiv:2108.07542 [physics.class-ph] .
  * Gutzwiller (1990) M. C. Gutzwiller,  _Chaos in Classical and Quantum Mechanics_ (Springer-Verlag, New York, NY, 1990).
  * Gray, Karl, and Novikov (1996) C. G. Gray, G. Karl, and V. A. Novikov, “Direct use of variational principles as an approximation technique in classical mechanics,” American Journal of Physics 64, 1177–1184 (1996).
  * Note (4) The black dotted curve in the figure is known as the caustic, and it is an envelope that divides spacetime into regions that are accessible from the given initial conditions, and regions that aren’t (see Ref. 30). Generally speaking, the caustic is composed of the kinetic foci associated with the family of trajectories that are dynamically allowed within a given potential (i.e. the family of solutions to the system’s equation or equations of motion). I will not be concerned with caustics in the rest of this paper, because the examples I consider in Secs. IV.1 and IV.2, namely harmonic oscillators with time-independent and time-dependent damping and frequency, do not exhibit caustics. Harmonic oscillators are special in that trajectories with arbitrarily large differences in their initial velocities reach the same kinetic focus, so for these systems the caustic curve collapses to a point Gray and Taylor (2007).
  * Morse, Marston (1934) Morse, Marston,  _The Calculus of Variations In The Large_ (The American Mathematical Society, New York, 1934).
  * Levit and Smilansky (1977) S. Levit and U. Smilansky, “A new approach to Gaussian path integrals and the evaluation of the semiclassical propagator,” Annals of Physics 103, 198–207 (1977).
  * Schumacher, Benjamin and Westmoreland, Michael D. (2010) Schumacher, Benjamin and Westmoreland, Michael D.,  _Quantum Processes, Systems, and Information_ (Cambridge University Press, New York, 2010).
  * Shankar, Ramamurti (1994) Shankar, Ramamurti,  _Principles of Quantum Mechanics_ (Springer, New York, 1994).
  * Note (5) This result was also obtained by G. Leitmann in Ref. 15, with an action of the form $S=\intop\ilimits@e^{2ct}\left[\dot{x}^{2}-\frac{1}{2}\omega_{0}^{2}x^{2}\right]dt$, where my $a$ equals his $2c$.
  * Lewis (1967) H. R. Lewis, “Classical and quantum systems with time-dependent harmonic-oscillator-type hamiltonians,” Physical Review Letters 18, 510 (1967).
  * Robles-Pérez (2017) S. Robles-Pérez, “Invariant vacuum,” Physics Letters B 774, 608–615 (2017), arXiv:1706.05474 [hep-th] .
  * Finelli, Vacca, and Venturi (1998) F. Finelli, G. P. Vacca, and G. Venturi, “Chaotic inflation from a scalar field in nonclassical states,” Physical Review D 58, 103514 (1998), arXiv:gr-qc/9712098 [gr-qc] .
  * Liu, Torres, and Wang (2018) Q. Liu, P. J. Torres, and C. Wang, “Contact Hamiltonian dynamics: Variational principles, invariants, completeness and periodic behavior,” Annals of Physics 395, 26–44 (2018).
  * Note (6) When I called this routine, I used the Sequential Least Squares Programming method. Documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html.
  * Note (7) https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html.
  * Note (8) Like the undamped harmonic oscillator, the TDHO is a special case in that trajectories having arbitrarily different initial velocities will intersect at the kinetic focus; it is not necessary for the difference between the initial velocities of two trajectories to be small.
  * Hunter (2007) J. D. Hunter, “Matplotlib: A 2D Graphics Environment,” Computing in Science and Engineering 9, 90–95 (2007).
  * Harris  _et al._ (2020) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” Nature 585, 357–362 (2020), arXiv:2006.10256 [cs.MS] .
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1. 0 Contributors, “SciPy 1.0: fundamental algorithms for scientific computing in Python,” Nature Methods 17, 261–272 (2020), arXiv:1907.10121 [cs.MS] .
  * Lazo and Krumreich (2014) M. J. Lazo and C. E. Krumreich, “The action principle for dissipative systems,” Journal of Mathematical Physics 55, 122902 (2014), arXiv:1412.5109 [math-ph] .
  * Tartaglia (1983) A. Tartaglia, “Non-conservative forces, lagrangians and quantisation,” European Journal of Physics 4, 231–234 (1983).
  * Kobe, Reali, and Sieniutycz (1986) D. H. Kobe, G. Reali, and S. Sieniutycz, “Lagrangians for dissipative systems,” American Journal of Physics 54, 997–999 (1986).
  * Ray (1979) J. R. Ray, “Lagrangians and systems they describe-how not to treat dissipation in quantum mechanics,” American Journal of Physics 47, 626–629 (1979).
  * Greenberger (1979) D. M. Greenberger, “A critique of the major approaches to damping in quantum theory,” Journal of Mathematical Physics 20, 762–770 (1979).
  * Edwards (1979) I. K. Edwards, “Quantization of inequivalent classical Hamiltonians,” American Journal of Physics 47, 153–155 (1979).
  * Cheng and Hong (2022) W. Cheng and J. Hong, “Local strict singular characteristics: Cauchy problem with smooth initial data,” Journal of Differential Equations 328, 326–353 (2022), arXiv:2103.06217 [math.AP] .