Levinson–Smith Dissipative Equations and Geometry of GENERIC Formalism and Contact Hamiltonian Mechanics

Abstract

We apply Jacobi’s Last Multiplier theory to construct the non-standard Lagrangian and Hamiltonian structures for the Levinson–Smith equations satisfying the Chiellini integrability condition. Then after a brief exposition of the contact geometry, we explore its connection with the non-standard Hamiltonian structures. We present the formulation of the Levinson–Smith equation in terms of General Equation for the Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) method and also study the gradient-type flow. We give a geometric formulation of GENERIC and apply this to general Levinson–Smith equations.

1 Introduction

The usual Liénard equation (Liénard 1928)

$$\begin{aligned} \ddot{x} +f(x){\dot{x}}+g(x)=0, \end{aligned}$$

(1)

for which a particular example, as indicated in Levinson and Smith (1942), is the Van der Pol equation (van der Pol 1926)

$$\begin{aligned} \ddot{x}+\mu (x^2-1) \dot{x}+x=0, \end{aligned}$$

(2)

is an autonomous second-order nonlinear differential equation playing an important rôle in modeling various phenomena in many different areas of the applied sciences. It has received a great deal of attention over the years (Harko et al. 2014; Bagchi et al. 2015; Kudryashov and Sinelshchikov 2015, 2016; Choudhury et al. 2017; Mitra et al. 2024) and has been generalized by Levinson and Smith (Levinson and Smith 1942) by considering a more general function $F(x,\dot{x})$ in the damping term:

$$\begin{aligned} \ddot{x}+F(x,\dot{x}) \dot{x}+g(x)=0, \end{aligned}$$

(3)

an equation similar to (1) but with f(x) replaced by $F(x, {\dot{x}})$. A particular case of the Levinson–Smith equation which is considered here is when $F(x, {\dot{x}})$ is an affine function of ${\dot{x}}$, i.e., $F(x, {\dot{x}})=p(x)\,\dot{x}+q(x)$, the equation being (Kudryashov and Sinelshchikov 2017)

$$\begin{aligned} \ddot{x} + p(x)\,\dot{x}^2+q(x){\dot{x}}+g(x)=0, \end{aligned}$$

(4)

which for $p\equiv 0$ reduces to the standard Liénard equation (Liénard 1928), while when $q\equiv 0$ it reduces to the quadratic Liénard equation (Kudryashov and Sinelshchikov 2015; Sabatini 2004; Tiwari et al. 2013). Moreover, it has been proved in Cariñena et al. (2023) that if a second-order differential equation is linearizable under a generalized Sundman transformations, then it should be of the form (4).

The integrability properties of a system of a system of differential equations, or of the corresponding vector field, are consequences of the existence of compatible geometric structure and the best knows case is complete integrability in the Arnold-Liouville sense, usually described in a symplectic o Poisson setting, but other possibilities are also relevant (Cariñena 2024). We are here intererested in the use of contact structures in the study of Levinson–Smith equation.

Contact geometry has been quite often used in physics, not only in optics but also in thermodynamics (Mrugala et al. 1991; Bravetti et al. 2015; Esen et al. 2022a, b). The geometry of equilibrium thermodynamics is extremely rich and based on contact geometry. In fact the phase space structure of the Liénard equation for a non-standard Hamiltonian framework bears a close resemblance with this geometry. In particular, the phase space turns out to be equipped with pseudo-Riemannian structures, whose restrictions to Legendre sub-manifolds define the thermodynamic equilibrium.

Unfortunately, the use of contact structures in geometric mechanics is less known at least till the last decade, and this suggests us to add in this paper a quick summary of the most relevant results of contact mechanics, both in the Hamiltonian and Lagrangian settings (Bravetti et al. 2017; de León et al. 2020).

Another important example of Eq. (4) appears when studying the dynamics of a 1-dimensional position-dependent mass motion under the action of a natural system, the equation including a ${\dot{x}}^{2}$-dependent additional term, which is a prototypical example of the Levinson–Smith equation (Mitra et al. 2024). Physically Eq. (4) may be regarded as a dissipative position-dependent mass system (Cariñena and Santos 2021). Such equations appear in many areas of physics including the Nosé-Hoover dynamics Ezra (2004).

In Kamke (1977) Kamke proposes the following non-autonomous differential equation:

$$\begin{aligned} \ddot{x} - \frac{f'(x)}{f(x)}{\dot{x}}^{2} + g(t){\dot{x}} + h(t)f(x) = 0, \end{aligned}$$

(5)

where f, g, h are arbitrary functions, and considers this equation as a particular instance of one of the form

$$\begin{aligned} \ddot{x} + \phi (t,x,{\dot{x}}){\dot{x}} = \gamma (t,x), \end{aligned}$$

(6)

where

$$\begin{aligned} \phi (t,x,{\dot{x}}) = -\left( \frac{f'(x)}{f(x)}{\dot{x}} - g(t)\right) , \quad \gamma (t,x) = -h(t)f(x). \end{aligned}$$

The term $\gamma (t,x)$ on th right-hand side in (6) can be seen as a disturbance function while $\phi (t,x,{\dot{x}})$ plays the rôle of a nonlinear friction. The particular case of (6)

$$\begin{aligned} \ddot{x} - \left( \dot{x}-\frac{x}{t}\right) ^af(t, x)=0, \end{aligned}$$

(7)

was studied in Kamke (1977) and the solution for $a=2$ and $f(t,x)=-(t+x)^{-1}$ was found. Other particular cases can be found in Acevedo et al. (2023). Recently Demina (Demina 2023) formulated the necessary and sufficient conditions of Liouvillian integrability for non-degenerate near infinity polynomial Levinson–Smith differential equations and found a number of novel Liouvillian integrable subfamilies.

In this paper at first we illuminate the relations among contact manifolds and non-standard Hamiltonian structures of the Levinson–Smith equations. In particular, we study the differential geometric formulation of the Liénard equation satisfying Chiellini integrability condition using non-standard Lagrangian and Hamiltonian structures. It is noteworthy to say that the non-standard Hamiltonian of the Liénard equation follows from the reduction of the Levinson–Smith equation formalism. We explored the non-standard Hamiltonian structures of the Liénard equation and their connections with the contact mechanics in our earlier papers (Cariñena and Guha 2019, 2024) and in this paper we study the geometrical structure the General Equation for the Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) method and its applications to study generalized Levinson–Smith equations from the perspective of this geometrical structure.

This paper is organized as follows. In Sect. 2, we review the fundamental properties of contact geometry and associated contact Hamiltonian mechanics, and furthermore Jacobi structure and Jacobi bracket are also introduced in this section, as well as Lagrangian formulation is also developed and in particular it is shown that the integral curves of the contact Hamiltonian vector field defined by the Lagrangian are solutions of the known as system of Herglotz equations. Two examples have been added to illustrate the theory Sect. 3 that is devoted to present, by making use of the Jacobi last multiplier theory, an original non-standard Lagrangian formulation of the Levinson–Smith equation using Chiellini integrability condition. We illustrate our construction with position-dependent mass systems and relativistic damped oscillators, and the Hamiltonian formulation of the Levinson–Smith equation using contact geometry is also studied in this Sect. 4 and our construction is illustrated with a damped system. Finally we develop the study of the GENERIC formulation for Levinson–Smith equation with generalized damping term using metriplectic or GENERIC method in Sect. 4, and the results of the paper are discussed in the last section.

2 Digression: Contact Hamiltonian Lagrangian Formulation of Mechanics

Let $(M,\eta )$ be a $2n + 1$-dimensional strict contact manifold M with the contact 1-form $\eta$ satisfying the maximal non-integrability condition $\eta \wedge (d\eta )^{\wedge n} \ne 0$ (see e.g., Arnold (1989)). There exists a unique vector field E, called the Reeb vector field, such that

$$\begin{aligned} i(E)d\eta = 0, \qquad i(E)\eta =\eta (E)= 1, \end{aligned}$$

(8)

and this implies that ${\mathcal {L}}_E\eta =i(E)d\eta +d(i(E)\eta )=0$.

In this section, we will follow the notation of Cariñena and Guha (2024). A vector field X on the strict contact manifold $(M,\eta )$ is an infinitesimal contact automorphism if, and only if, there exists a differentiable function $\rho$ such that the Lie derivative of $\eta$ with respect to X is ${\mathcal {L}}_X\eta = \rho \, \eta$, which is equivalent to say that the 2n-dimensional distribution $\ker \eta$ is preserved by the flow of the vector field X. Then the vector field X is said to be a contact vector field. For instance, the Reeb vector field E is a contact vector field (with $\rho =0$). The coefficient $\rho$ is $\rho =E(\eta (X))$, because a simple contraction with E of both members of $i(X)d\eta +d(\eta (X)))\rho \, \eta$ leads to $\rho =E(\eta (X))$.

The Darboux theorem for a $(2n+1)$-dimensional strict contact manifolds (see e.g., Grabowska and Grabowski (2022) and references therein) is analogous to the Darboux theorem for symplectic manifolds and establishes the existence, around each point on M, of a set of local coordinates $(q^i,p_i,z)$, with $i=1,\ldots ,n$, such that the contact 1-form is given by

$$\begin{aligned} \eta =dz+ p_i\,dq^i, \end{aligned}$$

(9)

where hereafter summation on repeated indices is understood, with the Reeb vector field being then given by

$$\begin{aligned} E = \frac{\partial }{\partial z}. \end{aligned}$$

One can check that the map $\flat : {{\mathfrak {X}}}(M) \rightarrow \Omega ^1(M)$ defined by

$$\begin{aligned} \flat (X) = -i(X)d\eta + \eta (X)\eta , \end{aligned}$$

(10)

is a $C^{\infty }(M)$-module isomorphism, and its inverse map will be denoted $\#=\flat ^{-1}: \Omega ^1(M)\rightarrow {{\mathfrak {X}}}(M)$. In particular, the image under $\flat$ of the Reeb vector field defined by (8) is $\flat (E)=\eta$, i.e., $\#(\eta )=E$. Similarly, the restriction of $\flat$ on the subset of horizontal vector fields coincides with that of the map $-\widehat{d\eta }$, where $\widehat{d\eta }: {{\mathfrak {X}}}(M) \rightarrow \Omega ^1(M)$ is defined as $\widehat{d\eta }(X)=i(X)d\eta$.

The important point is that given a function $F\in C^\infty (M)$, there is a uniquely defined vector field $X_F$ such that

$$\begin{aligned} i(X_F) d\eta =-dF+ E(F) \eta ,\qquad \eta (X_F)=F. \end{aligned}$$

(11)

It is important to remark that the Lie derivative of $\eta$ with respect to $X_F$ is

$$\begin{aligned} {\mathcal {L}}_{X_F}\eta = i(X_F)d\eta +d(i(X_F)\eta )=-dF+ E(F) \eta +dF= E(F)\eta , \end{aligned}$$

(12)

and hence $X_F$ is a contact vector field. Moreover, the first condition on (11) can be replaced by ${\mathcal {L}}_{X_F}\eta = E(F)\eta$, because then $i(X_F)d\eta ={\mathcal {L}}_{X_F}\eta -d(i(X_F)\eta =E(F)\eta -dF$, i.e., we recover the first condition.

Finally, as ${\mathcal {L}}_{X_F}(d\eta ) =d({\mathcal {L}}_{X_F}\eta ) =d(E(F)\eta )= E(F)\, d\eta + d(E(F))\wedge \eta$, we see that

$$\begin{aligned} {\mathcal {L}}_{X_F}\big (\eta \wedge (d\eta )^{\wedge n} \big ) = (n+ 1){E}(F) \big (\eta \wedge (d\eta )^{\wedge n} \big ). \end{aligned}$$

(13)

Note also that the image under $\flat$ of the contact vector field $X_F$ is given by

$$\begin{aligned} \flat (X_F) & =- i(X_F)d\eta + \eta (X_F)\eta =- i(X_F)d\eta + F\, \eta \Longleftrightarrow X_F\\ & =\#(- i(X_F)d\eta )+F\, {\mathcal {R}}_{\eta }, \end{aligned}$$

and rewriting the first condition in (11) as follows:

$$\begin{aligned} dF= - i(X_F)d\eta + {\mathcal {R}}_\eta (F)\,\eta , \end{aligned}$$

we see that

$$\begin{aligned} \#(dF)= X_F+({\mathcal {R}}_\eta (F)-F)\, {\mathcal {R}}_\eta . \end{aligned}$$

(14)

The vector field $X_F$ satisfying (11) is called the contact Hamiltonian vector field associated to the function F. In particular, when $F=\textbf{1}_M$ is the constant function equal to 1, we see that $X_{\textbf{1}_M}=E$. Moreover, if $F_1,F_2\in C^\infty (M)$, we see that $X_{F_1+F_2}=X_{F_1}+X_{F_2}$, while for each real constant c as $X_{cF}=cX_F$, hence we have $X_{c{\textbf{1}}_M}= c\,E$.

Given a function $F\in C^\infty (M)$, the coordinate expression of the vector field $X_F$ in a Darboux coordinate system $(q^i,p_i,z)$ is de León et al. (2020); Jin and Yan (2021); de León and Sardón (2017):

$$\begin{aligned} X_F=\frac{\partial F}{\partial p_i}\frac{\partial }{\partial q^i}+ \left( -\frac{\partial F}{\partial q^i}+p_i\frac{\partial F}{\partial z}\right) \frac{\partial }{\partial p_i} +\left( -p_j\frac{\partial F}{\partial p_j}+F\right) \frac{\partial }{\partial z}. \end{aligned}$$

(15)

Its integral curves are the solutions of the so-called dissipation Hamilton equations (Bravetti et al. 2017; de León and Sardón 2017):

$$\begin{aligned} \left\{ \begin{array}{rcl} {\displaystyle \frac{dq^i}{dt}} & =& {\displaystyle \frac{\partial F}{\partial p_i}}\\ {\displaystyle \frac{dp_i}{dt}} & =& -{\displaystyle \frac{\partial F}{\partial q^i}+p_i\frac{\partial F}{\partial z}}\\ {\displaystyle \frac{dz}{dt} }& =& F-p_j {\displaystyle \frac{\partial F}{\partial p_j}}\end{array} \right. \qquad . \end{aligned}$$

(16)

Furthermore, a contraction with the vector field $X_F$ of the first relation in (11) shows that $X_F(F)=F\, E(F)$, and then, $X_F(F^2)=2\,F^2\, E(F)$, and more generally, one can check by induction that, for any natural number k, $X_F(F^k)=k\,F^k\, E(F)$.

On the other side, from $d(F_1\,F_2)=F_1\,dF_2+F_2\,dF_1$ and $E(F_1\,F_2)=F_1\,E(F_2)+F_2\,E(F_1)$, we obtain

$$\begin{aligned} d(F_1\, F_2)=F_1(-i(X_{F_2})d\eta +E(F_2)\,\eta )+F_2(-i(X_{F_1})d\eta +E(F_1)\,\eta ) , \end{aligned}$$

and then,

$$\begin{aligned} d(F_1\, F_2)=-i(F_1\,X_{F_2}+F_2\,X_{F_1})d\eta +E(F_1\,F_2)\,\eta , \end{aligned}$$

and using definition (11) for $X_{F_1\,F_2}$, we see that there exists a real coefficient $\alpha$ such that

$$\begin{aligned} X_{F_1\,F_2}= F_1\, X_{F_2}+ F_2X_1+\alpha \, E. \end{aligned}$$

Such coefficient $\alpha$ is fixed by the condition $\eta (X_{F_1\,F_2})=F_1\,F_2$ as $\alpha =-F_1\,F_2$, and therefore,

$$\begin{aligned} X_{ F_2\, F_1} = F_2\,X_{ F_1}+ F_1\,X_{ F_2}- F_1\, F_2\,E. \end{aligned}$$

(17)

In the particular case $F_1=F_2$ in (17), we obtain $X_{F^2}=2\,F\, X_F-F^2\,E$.

2.1 Jacobi Structure and Jacobi Bracket

A strictly contact structure $(M,\eta )$ defines a Jacobi structure on M. We recall that a Jacobi structure on M is a pair $(\Lambda ,E)$ consisting of a vector field $E \in {{\mathfrak {X}}}(M)$ and a skew-symmetric bivector (called Jacobi bivector) field $\Lambda$ on M satisfying

$$\begin{aligned} [E,\Lambda ] =0, \qquad \hbox { and } \qquad [\Lambda ,\Lambda ] = 2E \wedge \Lambda , \end{aligned}$$

(18)

where $[\cdot , \cdot ]$ stands for the Schouten bracket (Ibort et al. 1997; de León et al. 2001). If $E= 0$, then $[\Lambda ,\Lambda ] = 0$, and thus we obtain as a particular case that of a Poisson bivector field giving rise to a Poisson structure.

A Jacobi structure $(\Lambda ,E)$ on M has associated a Jacobi bracket (see e.g., de León and Lainz (2020, 2019)) given by

$$\begin{aligned} \{F_2,F_1\}=\Lambda (dF_2,dF_1)+F_2\,E(F_1)-F_1\,E(F_2). \end{aligned}$$

The Jacobi structure determined by a contact manifold $(M,\eta )$ is defined as follows de León et al. (2001). First, the bivector field $\Lambda$ is given by

$$\begin{aligned} \Lambda (\alpha ,\beta )= d\eta (\flat ^{-1}(\alpha ), \flat ^{-1}(\beta ))=d\eta (\#(\alpha ), \#(\beta )), \end{aligned}$$

(19)

where $\flat : {{\mathfrak {X}}}(M) \rightarrow \bigwedge ^1(M)$, is the above mentioned of $C^{\infty }(M)$-module isomorphism defined by (10). Then, as relation (14) establishes that $\#(dF)$ and $X_F$ differ in a vertical vector field, we can also rewrite, for $\alpha =dF_2$ and $\beta =dF_1$, $\Lambda (dF_2,dF_1)= d\eta (X_{F_2},X_{F_1})$.

The other ingredient in the definition of Jacobi structure, the vector field E, is just the Reeb vector field of $(M,\eta )$, which according to (8), turns out to be $E ={\mathcal {R}}_\eta = \flat ^{-1}(\eta )=\#(\eta )$.

Then, a ‘Jacobi bracket’ on a strictly contact manifold $(M,\eta )$ can be defined by:

$$\begin{aligned} \{F_2,F_1\}= \Lambda (dF_2,dF_1)+F_2\, {\mathcal {R}}_\eta (F_1)- F_1\, {\mathcal {R}}_\eta (F_2), \end{aligned}$$

(20)

which is obviously ${\mathbb {R}}$-bilinear and skew-symmetric, and using that

$$\begin{aligned} \Lambda (dF_2,dF_1)=d\eta (X_{F_2},X_{F_1}) =(i(X_{F_2})d\eta )(X_{F_1}) =(-dF_2+{\mathcal {R}}_\eta (F_2)\,\eta )(X_{F_1}), \end{aligned}$$

namely

$$\begin{aligned} \Lambda (dF_2,dF_1)=-X_{F_1}(F_2) + F_1\, {\mathcal {R}}_\eta (F_2), \end{aligned}$$

(21)

and according to (12), it can also be rewritten as

$$\begin{aligned} \{F_2,F_1\}=-i(X_{F_1})(dF_2)+ F_2\, {\mathcal {R}}_\eta (F_1), \end{aligned}$$

(22)

or equivalently, using skew-symmetry of the bracket,

$$\begin{aligned} \{F_2,F_1\}=i(X_{F_2})(dF_1)-F_1\, {\mathcal {R}}_\eta (F_2). \end{aligned}$$

(23)

In particular, if $\textbf{1}_M$ denotes the constant function equal to one,

$$\begin{aligned} \{F,\textbf{1}_M\}=- {\mathcal {R}}_\eta (F). \end{aligned}$$

(24)

The Jacobi bracket does not satisfy the Leibnitz rule. For instance, the preceding relation shows that

$$\begin{aligned} \{F_2,F_1\}=\{F_2,\textbf{1}_M\,F_1\} \ne \textbf{1}_M\,\{F_2,F_1\}+F_1\,\{F_2,\textbf{1}_M\}. \end{aligned}$$

More explicitly, (22) means that

$$\begin{aligned} \{\cdot ,F\}=-X_F(\cdot )+{\mathcal {R}}_\eta (F)\, \cdot , \end{aligned}$$

(25)

i.e.,

$$\begin{aligned} \{G, F\}=-X_F(G)+G\, {\mathcal {R}}_\eta (F), \quad \forall G\in C^\infty (M), \end{aligned}$$

(26)

and then $\{F_3,F_2\,F_1\}\ne \{F_3,F_2\}\,F_1+F_2\, \{F_3,F_1\}$, because of the last term of the right-hand side in the preceding relation, but

$$\begin{aligned} \{F_3,F_2\,F_1\}=\{F_3,F_2\}\,F_1+F_2\, \{F_3,F_1\}+F_1\,F_2\,{\mathcal {R}}_\eta (F_3), \end{aligned}$$

i.e.,

$$\begin{aligned} \{F_3,F_2\,F_1\}=\{F_3,F_2\}\,F_1+F_2\, \{F_3,F_1\}-F_1\,F_2\,\{F_3,\textbf{1}_M\}. \end{aligned}$$

(27)

Proposition 1

The commutator $[X_{F_2},X_{F_1}]$ of two contact Hamiltonian vector fields is a contact Hamiltonian vector field, and more explicitly, for each pair of functions, $[X_{F_2},X_{F_1}]=X_{\{F_2,F_1\}}$.

It is clear that if $X_{F_2}$ and $X_{F_1}$ preserve the contact distribution, its commutator $[X_{F_2},X_{F_1}]$ also preserves such contact distribution and then it is a contact vector field. In order to determine the function defining such contact Hamiltonian vector field, we can make use of the relation

$$\begin{aligned} ({\mathcal {L}}_{X_{F_2}}\circ i(X_{F_1})- i(X_{F_1})\circ {\mathcal {L}}_{X_{F_2}})\eta =i([X_{F_2},X_{F_1}])\eta , \end{aligned}$$

and recalling that ${\mathcal {L}}_{X_{F_2}}\eta ={\mathcal {R}}_\eta (F_2)\, \eta$, we obtain

$$\begin{aligned} X_{F_2}(F_1)-F_1\,{\mathcal {R}}_ \eta (F_2)=i([X_{F_2},X_{F_1}])\eta , \end{aligned}$$

which, using (23), as $X_{F_2}(F_1)-{\mathcal {R}}_\eta (F_2)\,F_1=\{F_2,F_1\}$, shows that $\eta ([X_{F_2},X_{F_1}])=\{F_2,F_1\}$ and hence the contact Hamiltonian vector field corresponding to $[X_{F_2},X_{F_1}]$ is

$$\begin{aligned} [X_{F_2},X_{F_1}]=X_{\{F_2,F_1\}}. \end{aligned}$$

(28)

$\square$

This relation can be used to prove that the composition law $\{\cdot ,\cdot \}$ satisfies Jacobi identity. In fact, we recall that the contact vector field $X_F$ is zero if, and only if, $F=0$, and we can compute, for any three functions $F_i\in C^\infty (M)$, $i=1,2,3$,

$$\begin{aligned} X_{\{F_i,\{F_j,F_k\}\}}=\left[ X_{F_i},X_{\{F_j,F_k\}}\right] =\left[ X_{F_i},[X_{F_j},X_{F_k}]\right] , \end{aligned}$$

and then, Jacobi identity for $\{\cdot ,\cdot \}$ is consequence of Jacobi identity for the commutator of vector fields.

Therefore $\{\cdot ,\cdot \}$ endows the set $C^\infty (M)$ with a real Lie algebra structure, and (28) shows that the association of a function with its corresponding contact vector field is a homomorphism of Lie algebras.

2.2 Contact Hamiltonian

In a Darboux coordinate system the expressions of $\Lambda$ and $\{F_2,F_1\}$ are:

$$\begin{aligned} \begin{array}{rcl} \Lambda & =& {\displaystyle \left( \frac{\partial }{\partial q^i}\wedge \frac{\partial }{\partial p_i}+p_i\, \frac{\partial }{\partial z}\wedge \frac{\partial }{\partial p_i}\right) } \\ & & \\ \{F_2,F_1\}& =& -\left( {\displaystyle \frac{\partial F_2}{\partial q^i}\,\frac{\partial F_1}{\partial p_i}- \frac{\partial F_2}{\partial p_i}\,\frac{\partial F_1}{\partial q^i}+p_i\frac{\partial F_1}{\partial z}\frac{\partial F_2}{\partial p_i}- p_i\frac{\partial F_1}{\partial p_i}\frac{\partial F_2}{\partial z}}\right) \\ & & \quad -F_1{\displaystyle \frac{\partial F_2}{\partial z} +F_2\frac{\partial F_1}{\partial z}} \end{array} \end{aligned}$$

(29)

which exhibits that $X_F(G)= \{G,F\}+E(F)G$.

The contact (horizontal) distribution is defined by

$$\begin{aligned} {{{\mathcal {H}}}} = \hbox { ker}(\eta ) = \{{\xi } \in TM \mid \eta ({\xi }) = 0 \}, \end{aligned}$$

(30)

and the line bundle $V_{E}$ or vertical distribution is generated by the Reeb vector field E:

$$\begin{aligned} V_{E} = \hbox { span } ({E})= \hbox { ker}(d\eta ), \quad \hbox { where } \,\,\,\,\, {E} = \frac{\partial }{\partial z}. \end{aligned}$$

(31)

Thus the tangent bundle of M has the natural decomposition, $TM = V_{E} \oplus {{{\mathcal {H}}}}$. The subbundle of rank 2n ${{{\mathcal {H}}}}$ $(q^i,p_i,z)$, is generated in a Darboux coordinate system by the set of vector fields

$$\begin{aligned} \left\{ X_{q^i}=-\frac{\partial }{\partial p_i},\, X_{p_i}=-p_i\frac{\partial }{\partial z} + \frac{\partial }{\partial q^i}\ \Big |\ i=1,\ldots ,n \right\} , \end{aligned}$$

satisfying the relations, for $i,j=1,\ldots ,n$,

$$\begin{aligned} & [X_{q^i}, X_{q^j}] =0,\quad [X_{q^i}, X_{p_j}] = \delta _{j}^{i}{E}, \quad [ X_{q^i}, {E}] = 0,\quad [X_{p_i}, X_{p_j}] =0,\\ & [ X_{p_i}, {E}] = 0. \end{aligned}$$

If in the expression (11) for the contact vector field $X_F$ we replace the generic function F by the standard Hamiltonian notation $H(q^i,p_i,z)$ we see that an integral curve $\big (q^i(t),p_i(t),z(t) \big )$ of $X_H$ is a solution of the set of dissipative Hamiltonian equations

$$\begin{aligned} \left\{ \begin{array}{rcl} {\displaystyle \frac{dq^i}{dt}} & =& {\displaystyle \frac{\partial H}{\partial p_i}}=-\{q^i,H\} + q^i\,{\displaystyle \frac{\partial H}{\partial z}}\\ {\displaystyle \frac{dp_i}{dt}} & =& -{\displaystyle \frac{\partial H}{\partial q^i}+p_i\frac{\partial H}{\partial z}}=-\{p_i,H\}+p_i {\displaystyle \frac{\partial H}{\partial z}}\\ {\displaystyle \frac{dz}{dt} }& =& H-p_j {\displaystyle \frac{\partial H}{\partial p_j}}=-\{z,H\}+ z\,{\displaystyle \frac{\partial H}{\partial z}}\end{array} \right. \qquad , \end{aligned}$$

(32)

which corresponds to (25): $X_H(\cdot )=-\{\cdot , H\}+ R_\eta (H)\cdot$.

These systems of equations can be regarded as instances of Hamiltonian systems on Jacobi manifolds. Unlike the Poisson (or conservative) case, the vector field $X_H$ satisfies $X_H(H) = H{E}(H)$. Recall also that we have proved at (12), with $F=H$, that the Lie derivative of $\eta$ with respect to $X_H$ is ${\mathcal {L}}_{X_H}\eta = {E}(H)\eta$, which yields to

$$\begin{aligned} {\mathcal {L}}_{X_H}\big (\eta \wedge (d\eta )^{\wedge n} \big ) = (n+ 1){E}(H) \big (\eta \wedge (d\eta )^{\wedge n} \big ). \end{aligned}$$

(33)

2.3 Lagrangian Formulation of Contact Mechanics

The Lagrangian formulation of contact mechanics makes use of the Herglotz variational principle to study dissipative systems with a new variable s (de León and Lainz 2020, 2019). Consider a function $L\in C^\infty (TQ \times \mathbb {R})$ that we will still call Lagrangian. Its expression in local coordinates $(q^i,v^i,s)$ induced from local coordinates $(q^i)$ on Q, $(q^i,v^i)$ being the corresponding coordinates on TQ and s the natural coordinate on ${\mathbb {R}}$, L is of the form $L(q^i,v^i,s)$, where $L: TQ \times \mathbb {R} \rightarrow \mathbb {R}$.

Recall that the usual tangent structure of TQ (see e.g., Cariñena et al. (2023) and references therein) is characterized by the Liouville dilation vector field $\Delta$ and the vertical endomorphism S. In local coordinates on TQ induced from local coordinates in Q their local expressions are:

$$\begin{aligned} \Delta =v^i\, \frac{\partial }{\partial v^i},\qquad S= \frac{\partial }{\partial v^i}\otimes d q^i. \end{aligned}$$

They can be extended to $TQ \times \mathbb {R}$ with trivial action on basic functions (functions on $\mathbb {R}$) and vector fields on $TQ\times \mathbb {R}$ that are projectable on $\mathbb {R}$, and these extensions will be denoted by the same symbols $\Delta$ and S.

Then, given such a Lagrangian function we define the 1-form $\theta _L=dL\circ S$ on $TQ \times \mathbb {R}$ as the (pull-back on $TQ \times \mathbb {R}$ of the corresponding) Lagrangian 1-form on TQ, that in the above-mentioned local coordinates is written as

$$\begin{aligned} \theta _L= \frac{\partial L}{\partial v^i}\,dq^i. \end{aligned}$$

Note that the 1-form $\eta _L$ on $TQ \times \mathbb {R}$ defined by

$$\begin{aligned} \eta _L = ds + \theta _L = ds + \frac{\partial L}{\partial v^i}\,dq^i, \end{aligned}$$

(34)

is such that, as

$$\begin{aligned} d\eta _L & = d\theta _L = - dq^i \wedge d\left( \frac{\partial L}{\partial v^i} \right) = - \frac{\partial ^2 L}{\partial q^j \partial v^i}dq^i \wedge dq^j - \frac{\partial ^2 L}{\partial v^j \partial v^i}dq^i \wedge dv^j \\ & \quad \ - \frac{\partial ^2 L}{\partial s \partial v^i}dq^i \wedge ds, \end{aligned}$$

the 1-form $\eta _L$ satisfies the maximal non-integrability condition

$$\begin{aligned} \eta _L \wedge (d\eta _L)^{\wedge n} \ne 0, \end{aligned}$$

(35)

if, and only if, $\det W\ne 0$, where W denotes the Hessian matrix with elements

$$\begin{aligned} W_{ij}=\dfrac{\partial ^2L}{\partial v^i\,\partial v^j}, \end{aligned}$$

and this regularity condition for L, which. will be assumed hereafter, is the condition in terms of L for the 1-form $\eta _L$ to be a co ntact 1-form on $TQ \times \mathbb {R}$.

Proposition 2

The coordinate expression of the Reeb vector field $R_L$ is given by

$$\begin{aligned} {R_L} = \frac{\partial }{\partial s} - W^{ij}\frac{\partial ^2 L}{\partial v^i \partial s} \frac{\partial }{\partial v^j}, \end{aligned}$$

(36)

where $W^{ij}$ are the elements of the inverse of the Wronskian matrix, i.e., $W^{ij}W_{jk}=\delta ^i_k$.

Proof

In fact, if the local expression of the Reeb vector field $R_L$ is

$$\begin{aligned} R_L=E^k\frac{\partial }{\partial q^k}+F^k\frac{\partial }{\partial v^k}+E^s\frac{\partial }{\partial s}, \end{aligned}$$

then

$$\begin{aligned} i(R_L)d\eta _L & =-E^k\left( \left( \frac{\partial ^2L}{\partial q^k\partial v^j}-\frac{\partial ^2L}{\partial q^j\partial v^k}\right) dq^j+\frac{\partial ^2L}{\partial v^k\partial v^j} dv^j\right) +F^k\frac{\partial ^2L}{\partial v^k\partial v^j} dq^j \\ & \quad \ + E^s \frac{\partial ^2 L}{\partial s\partial v^i}dq^i, \end{aligned}$$

and therefore, because of the assumed regularity of L, the condition $i(R_L)d\eta _L= 0$ implies that $E^k=0$, and furthermore,

$$\begin{aligned} W_{ij} F^j+E^s\frac{\partial ^2 L}{\partial v^i \partial s}=0,\quad j=1,\ldots ,n, \end{aligned}$$

that, together with the condition $\eta _L(R_L)=1$, imply the mentioned form of the Reeb vector field $R_L$. $\square$

Proposition 3

If the Lagrangian function $L\in C^\infty (TQ\times {\mathbb {R}})$ is regular, then the coordinate expression of the contact Hamiltonian vector field defined by the energy function of the system, $E_L$, defined on the strictly contact manifold $(TQ\times {\mathbb {R}},\eta _L)$ by

$$\begin{aligned} E_L=\Delta (L)-L=v^i\, \frac{\partial }{\partial v^i}-L, \end{aligned}$$

(37)

is given by

$$\begin{aligned} \Gamma _L=X_{E_L}=v^i\frac{\partial }{\partial q^i}+B^i\frac{\partial }{\partial v^i}-L\frac{\partial }{\partial s}, \end{aligned}$$

(38)

where the coefficients $B^i$ satisfy the following conditions:

$$\begin{aligned} B^k\frac{\partial ^2L}{\partial v^k\partial v^j} -L \frac{\partial ^2 L}{\partial s\partial v^j}+v^k\frac{\partial ^2L}{\partial q^k\partial v^j}-\frac{\partial L}{\partial q^j}-\frac{\partial L}{\partial s}\frac{\partial L}{\partial v^j}=0, \quad j=1,\ldots ,n. \end{aligned}$$

Proof

As the energy $E_L$ of the system is defined by (37), when L is regular, the dynamics in the Lagrangian formalism is given by the contact Hamiltonian vector field defined by $E_L$, i.e., it is the univocally defined vector field $\Gamma _L=X_{E_L}$ such that

$$\begin{aligned} \eta _L(\Gamma _L)= E_L,\quad i(\Gamma _L)d\eta _L=-dE_L+R_L(E_L)\,\eta _L. \end{aligned}$$

(39)

Note that if

$$\begin{aligned} \Gamma _L=A^i\frac{\partial }{\partial q^i}+B^i\frac{\partial }{\partial v^i}+C\frac{\partial }{\partial s}, \end{aligned}$$

(40)

the first condition of (39) implies that, as

$$\begin{aligned} \eta (\Gamma _L)=\left( ds+\frac{\partial L}{\partial v^k}\, dq^k\right) \left( A^i\frac{\partial }{\partial q^i}+B^i\frac{\partial }{\partial v^i}+C\frac{\partial }{\partial s}\right) =C+A^i\frac{\partial L}{\partial v^i}= v^i\frac{\partial L}{\partial v^i}-L, \nonumber \\ \end{aligned}$$

(41)

we must have

$$\begin{aligned} C+A^i\frac{\partial L}{\partial v^i}= v^i\frac{\partial L}{\partial v^i}-L, \end{aligned}$$

while for $i(\Gamma _L)d\eta _L=i(\Gamma _L)d\theta _L$ we have

$$\begin{aligned} i(\Gamma _L)d\eta _L & =-A^k\left( \left( -\frac{\partial ^2L}{\partial q^k\partial v^j}+\frac{\partial ^2L}{\partial q^j\partial v^k}\right) dq^j+\frac{\partial ^2L}{\partial v^k\partial v^j} dv^j\right) +B^k\frac{\partial ^2L}{\partial v^k\partial v^j} dq^j\nonumber \\ & +C \frac{\partial ^2 L}{\partial s\partial v^i}dq^i. \end{aligned}$$

(42)

On the other side, using the expression (36) for the Reeb vector field, we see that

$$\begin{aligned} R_L(E_L)=\frac{\partial }{\partial s}\left( v^i\, \frac{\partial L}{\partial v^i}-L\right) -W^{ij}\frac{\partial ^2L}{\partial s\partial v^i}\frac{\partial }{\partial v^j}\left( v^k\, \frac{\partial L}{\partial v^k}-L\right) , \end{aligned}$$

and from

$$\begin{aligned} \frac{\partial }{\partial v^j}\left( v^k\, \frac{\partial L}{\partial v^k}-L\right) =W_{jk}v^k, \quad j=1,\ldots ,n, \end{aligned}$$

we find that

$$\begin{aligned} R_L(E_L)=-\frac{\partial L}{\partial s}, \end{aligned}$$

and taking into account that

$$\begin{aligned} dE_L & =\frac{\partial L}{\partial v^j} dv^j+ v^k\frac{\partial ^2L}{\partial s\partial v^k}ds+ v^k\frac{\partial ^2L}{\partial q^j\partial v^k}dq^j+ v^k\frac{\partial ^2L}{\partial v^j\partial v^k}dv^j- \frac{\partial L}{\partial s} ds\\ & \quad \ -\frac{\partial L}{\partial q^j}dq^j- \frac{\partial L}{\partial v^j}dv^j, \end{aligned}$$

we see that

$$\begin{aligned} \begin{array}{rcl}dE_L-R_L(E_L){\displaystyle \left( ds+\frac{\partial L}{\partial v^k}\, dq^k\right) }& = {\displaystyle \left( v^k\frac{\partial ^2L}{\partial s\partial v^k}+\frac{\partial L}{\partial s}\right) ds +\left( v^j\frac{\partial ^2L}{\partial q^k\partial v^j}-\frac{\partial L}{\partial q^k}+\frac{\partial L}{\partial s}\frac{\partial L}{\partial v^k}\right) dq^k}\\ & + {\displaystyle v^k\frac{\partial ^2L}{\partial v^j\partial v^k}dv^j}\end{array} \end{aligned}$$

and comparing this expression with (42) we obtain that, for $j=1,\ldots ,n$,

$$\begin{aligned} \begin{array}{rl} & A^k{\displaystyle \left( -\frac{\partial ^2L}{\partial q^k\partial v^j}+\frac{\partial ^2L}{\partial q^j\partial v^k}\right) -B^k\frac{\partial ^2L}{\partial v^k\partial v^j} -C \frac{\partial ^2 L}{\partial s\partial v^j}} ={\displaystyle v^k\frac{\partial ^2L}{\partial q^j\partial v^k}-\frac{\partial L}{\partial q^j}+\frac{\partial L}{\partial s}\frac{\partial L}{\partial v^j}}\, \\ & A^k{\displaystyle \frac{\partial ^2L}{\partial v^k\partial v^j} = v^k\frac{\partial ^2L}{\partial v^j\partial v^k}}\\ & \end{array} \end{aligned}$$

(43)

and as L is assumed to be regular, the last equations show that $A^i=v^i$. In this case, the condition $\eta _L(\Gamma _L)=0$ gives for C the value $C=-L$. Replacing these values for $A^i$ and C on the first set of equations (43) we obtain the condition for the coefficients $B^i$:

$$\begin{aligned} B^k\frac{\partial ^2L}{\partial v^k\partial v^j} -L \frac{\partial ^2 L}{\partial s\partial v^j}+v^k\frac{\partial ^2L}{\partial q^k\partial v^j}-\frac{\partial L}{\partial q^j}-\frac{\partial L}{\partial s}\frac{\partial L}{\partial v^j}=0, \quad j=1,\ldots ,n. \end{aligned}$$

Therefore, the integral curves of $\Gamma _L$ are solutions of the known as system of Herglotz equations:

$$\begin{aligned} \left\{ \begin{array}{rl} & \dot{s}=-L\\ & {\displaystyle \frac{d}{dt}\left( \frac{\partial L}{\partial \dot{q}^i}\right) -\frac{\partial L}{\partial q^i}}={\displaystyle \frac{\partial L}{\partial \dot{q}^i}\frac{\partial L}{\partial s}} \end{array}\right. \qquad . \end{aligned}$$

(44)

$\Box$

Finally, remark that

$$\begin{aligned} i(\Gamma _L)\eta _L=(ds+\theta _L)(\Gamma _L)=-L+\theta _L(\Gamma _L)=-L+v^k\frac{\partial l}{\partial v}^k=-E_L, \end{aligned}$$

and then,

$$\begin{aligned} \flat _L(\Gamma _L)=-i(\Gamma _L)d\eta _L+\eta (\Gamma _L) =dE_L-R(E_L)\eta _L-E_L\eta _L=dE_L-(R(E_L)+E_L)\eta _L. \end{aligned}$$

We can also extend the Legendre transformation and define a map ${{{\mathcal {F}}}}L: TQ \times \mathbb {R} \rightarrow T^{*}Q \times \mathbb {R}$ that in local coordinates is given by

$$\begin{aligned} {{{\mathcal {F}}}}L(q^i, v^i,s) = \left( q^i,\frac{\partial L}{\partial v^i},s\right) , \end{aligned}$$

(45)

and this map relates the contact 1-form $\eta$ on $T^{*}Q \times \mathbb {R}$ that is the natural prolongation of the Liouville 1-form on $T^{*}Q \times \mathbb {R}$, to appear in the Hamiltonian framework, to the one in the Lagrangian framework,

$$\begin{aligned} ({{{\mathcal {F}}}}L)^{*}\eta = \eta _L. \end{aligned}$$

(46)

We see in (45) that the Lagrangian L is regular when ${{{\mathcal {F}}}}L$ is a (local) diffeomorphism. In this case $(TQ \times \mathbb {R}, \eta _L)$ is a strict contact manifold and if ${{{\mathcal {F}}}}L$ is a diffeomorphism we can use the inverse map of ${{{\mathcal {F}}}}L$ to define the Hamiltonian function, $H = E_L \circ ({{{\mathcal {F}}}}L)^{-1}\in C^\infty (T^{*}Q \times \mathbb {R})$, and the contact Hamiltonian vector field $X_H$ on the strict contact manifold $(T^{*}Q \times \mathbb {R},\eta )$ is ${{{\mathcal {F}}}}L$-related to the contact Hamiltonian vector field $X_{E_L}$ on $TQ \times \mathbb {R}$, namely $({{{\mathcal {F}}}}L)_{*}(X_{E_L}) = X_H$.

In the geometrical approach to mechanics, symplectic and related (presymplectic, Poisson, etc.) structures play a relevant rôle. However Riemannian geometry is also relevant in the geometric reformulation of classical systems of mechanical type, a Riemannian metric allowing us a definition, in an intrinsic way, of a kinetic energy function. Then, if (M, g) is a Riemannian manifold, i.e., g is a non-degenerate symmetric two-times covariant tensor field on M, we can define the kinetic energy $T_g\in C^\infty (TM)$ by

$$\begin{aligned} T_g(v)=\frac{1}{2}g(v,v),\qquad v\in TM, \end{aligned}$$

that can be considered as a Lagrangian for a free particle on M. Moreover if, additionally, a potential function $V\in C^\infty (M)$ is given, we can define a natural (standard or of mechanical type) Lagrangian $L\in C^\infty (TM)$ as follows (see e.g., Cariñena et al. (2012, 2014, 2022); Cariñena and Muñoz-Lecanda (2023)):

$$\begin{aligned} L_{g,V}(q,v)=T_g(v)-V(q)=\frac{1}{2} \,g_q(v,v)-(\tau _M^*V)(q,v)=\frac{1}{2} \,g_q(v,v)-V(q),\nonumber \\ \end{aligned}$$

(47)

where $T_g$ is the above kinetic energy. The simplest example is the 1-dimensional case which corresponds to the usually called position-dependent mass systems that will be used next to illustrate the general theory.

2.4 Position-Dependent Mass Dissipative Systems

The study of one-dimensional systems with position-dependent mass has received a lot of attention in recent times, with the non-commutativity of mass and momentum being a difficulty that needs to be taken into account in the quantization process (von Roos 1983; Lévy-Leblond 1995; Koç and Koca 2003; Cariñena et al. 2004, 2017). From the classical point of view the situation corresponds to systems for which the kinetic energy is defined from a non-Euclidean metric.

The kinetic energy function takes the form $T(x,v) = \dfrac{1}{2}m(x)\, v^{2}$, where m(x) is a positive differentiable function, and a standard Lagrangian L is of the form

$$\begin{aligned} L(x,v) = \frac{1}{2} \, m(x) v^2 - V(x). \end{aligned}$$

(48)

The Legendre transformation is given by ${{{\mathcal {F}}}}L(x,v)=(x,m(x)\,v)$, and the corresponding Hamiltonian of the position-dependent mass system is

$$\begin{aligned} H(x,p) = \frac{p^2}{2m(x)} + V(x). \end{aligned}$$

(49)

The Euler-Lagrange equation of motion for such a Lagrangian function (48) may be written as

$$\begin{aligned} m(x)\,\ddot{x} + \frac{1}{2} m^{\prime }(x)\,{\dot{x}}^2 + \frac{dV(x)}{dx} = 0. \end{aligned}$$

(50)

As an interesting instance, mass function $m(x) = \dfrac{1}{1 + \lambda x^2}$, $\lambda \in {\mathbb {R}}$, without singularities if $\lambda >0$, was proposed by Mathews and Lakshmanan (1974) in relation to the theory of Lagrangians of polynomial type in field theory and its quantum version was studied in Cariñena et al. (2004), while the generalization for multidimensional case was carried out in Cariñena et al. (2004) for the classical case and in Cariñena et al. (2007) for the quantum one.

We can include a Rayleigh damping term $\gamma \,{\dot{x}}$ in the preceding equation (Strutt 1873; López-Gordón 2021) and we obtain the position-dependent mass damped system

$$\begin{aligned} m(x)\,\ddot{x} + \frac{1}{2} m^{\prime }(x){\dot{x}}^2 + \gamma \,{\dot{x}} + \frac{dV(x)}{dx} = 0,\qquad \gamma \in {\mathbb {R}}, \end{aligned}$$

(51)

where $\gamma$ is the damping coefficient.

It is easy to check that the position-dependent mass damped system (51) can be constructed from contact Hamiltonian mechanics using as contact Hamiltonian the function

$$\begin{aligned} H_{\textrm{cpdmh}} = \frac{p^2}{2m(x)} + V(x) - \gamma \, s,\qquad \gamma \in {\mathbb {R}}. \end{aligned}$$

(52)

In fact, putting this explicit value of H in (32) we obtain as a first condition $\dot{x}=\dfrac{p}{m}$, and in the second one,

$$\begin{aligned} \dot{p}=-V'(x)+m'(x)\dot{x}^2+\frac{p^2}{2m^2}m'(x)-\gamma p, \end{aligned}$$

and therefore,

$$\begin{aligned} \frac{d}{dt}(m\,\dot{x})= m\ddot{x}+m'(x)\dot{x}^2= -V'(x)+\frac{1}{2} m'(x)\dot{x}^2-m\gamma \,\dot{x}, \end{aligned}$$

which coincides with (51).

2.5 Relativistic Damped Oscillator Equation

The simplest straightforward generalization of the damped oscillator equation of standard non-relativistic regime is to invoke the relativistic kinetic energy term. Recall that the Lagrangian of the one-dimensional free relativistic particle is defined as

$$\begin{aligned} L _0(\dot{x})= - m_0c^2\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}}, \end{aligned}$$

(53)

where $|\dot{x}|<c$, and then, as the variable x is cyclic, the relativistic momentum given by

$$\begin{aligned} p= \frac{\partial L_0}{\partial {\dot{x}}} = \frac{m_0{\dot{x}}}{\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}}}, \end{aligned}$$

(54)

is a constant of motion.

Observe also that the energy of the free relativistic particle is given by

$$\begin{aligned} E_{L_0}=\dot{x}\frac{m_0{\dot{x}}}{\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}}}+ m_0c^2\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}}=\frac{m_0c^2}{\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}}}. \end{aligned}$$

(55)

In the case of a one-dimensional relativistic harmonic oscillator, whose Lagrangian is

$$\begin{aligned} L (x,\dot{x})= - m_0c^2\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}} - \frac{1}{2}kx^2,\quad k>0, \end{aligned}$$

(56)

where $|\dot{x}|<c$, the equation of motion is

$$\begin{aligned} \frac{d}{dt}\left( \frac{m_0{\dot{x}}}{\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}}}\right) + k\,x = 0. \end{aligned}$$

(57)

According to Hutten’s description (Hutten 1965) the frequency decreases with the energy. After expanding the derivative term

$$\begin{aligned} m_0\ddot{x}\left( 1- \frac{{\dot{x}}^{2}}{c^2}\right) +\frac{m_0}{c^2}\dot{x}^2\ddot{x}+k\,x\left( 1- \frac{{\dot{x}}^{2}}{c^2}\right) ^{3/2}= m_0\ddot{x}+k\,x\left( 1- \frac{{\dot{x}}^{2}}{c^2}\right) ^{3/2}=0, \end{aligned}$$

we can approximate such equation by keeping the relativistic corrections up to order $1/c^2$, i.e.,

$$\begin{aligned} \left( 1- \frac{{\dot{x}}^{2}}{c^2}\right) ^{3/2}\approx 1-\frac{3}{2} \frac{{\dot{x}}^{2}}{c^2}, \end{aligned}$$

by

$$\begin{aligned} \ddot{x} + \omega ^2 x - \mu \,x\, {\dot{x}}^{2} = 0, \quad \textrm{with}\quad \omega ^2 = \frac{k}{m_0}, \quad \mu =\frac{3 \omega ^2}{2c^2}, \end{aligned}$$

(58)

which is a particular instance of (4).

Had we introduced Rayleigh friction term $\kappa \, {\dot{x}}$ in this framework (Strutt 1873), this would lead to the approximating damped relativistic oscillator equation

$$\begin{aligned} \ddot{x} - \mu \, x\, {\dot{x}}^{2} + \kappa \, {\dot{x}} + \omega ^2 x = 0, \end{aligned}$$

(59)

a particular instance of those whose linearization was studied in Sinelshchikov (2020).

Recall that as the Legendre transformation for the one-dimensional free particle is given by

$$\begin{aligned} {{{\mathcal {F}}}}L_0(x,v)=\left( x,\frac{\partial L_0}{\partial v}\right) = \left( x, \frac{m_0{\dot{x}}}{\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}}}\right) , \end{aligned}$$

for $|\dot{x}|<c$, such transformation is invertible, because from

$$\begin{aligned} p=\frac{m_0{\dot{x}}}{\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}}} \end{aligned}$$

we obtain

$$\begin{aligned} p^2+m_0^2c^2=m_0c^2\left( 1- \frac{{\dot{x}}^{2}}{c^2}\right) ^{-1} \Longleftrightarrow 1- \frac{{\dot{x}}^{2}}{c^2}=\frac{m_0^2c^2}{p^2+m_0^2c^2}, \end{aligned}$$

and hence the relation

$$\begin{aligned} \dot{x}=\frac{p}{m_0}\sqrt{1 - \frac{{\dot{x}}^{2}}{c^2}}=\frac{pc}{\sqrt{p^2+m_0^2c^2}}.\, \end{aligned}$$

Therefore the Hamiltonian for a free particle is

$$\begin{aligned} H_0(x,p)=p\,\frac{pc}{\sqrt{p^2+m_0^2c^2}} +m_0c^2\frac{m_0c}{\sqrt{p^2+m_0^2c^2}}=c\sqrt{p^2+m_0c^2}. \end{aligned}$$

We can also consider the particular case of a relativistic particle subject to a harmonic potential whose Hamiltonian is given by

$$\begin{aligned} H = c\sqrt{p^2 +m_0^2c^2} + \frac{1}{2}k\,x^2, \end{aligned}$$

(60)

and generalize this to the relativistic evolution of a particle with mass $m_0$ subject to a general potential V(x), and then.

$$\begin{aligned} H = c\sqrt{p^2 +m_0^2c^2} + V(x). \end{aligned}$$

(61)

The contact formalism for the standard free relativistic or for a particle under the harmonic oscillator potential, for instance, can be developed with a contact Hamiltonian defined in such a way that when the new variable s vanishes the contact Hamiltonians reduce to the corresponding standard Hamiltonians. This is not a Lorentz invariant expression, and then it is only relativistic in partial sense. This requires the extension to general relativity.

Consider, for instance, the relativistic ‘ansatz’

$$\begin{aligned} H_{c} = c\sqrt{p^2 +m_0^2c^2} + \frac{1}{2}k\,x^2- \gamma \, s. \end{aligned}$$

(62)

Then, the contact Hamiltonian equation yields

$$\begin{aligned} \left\{ \begin{array}{rcl} {\dot{x}} & =& {\displaystyle \frac{\partial H_c}{\partial p} = \frac{cp}{\sqrt{p^2 +m_0^2c^2}}}\\ {\dot{p}} & =& - {\displaystyle \frac{\partial H_c}{\partial x} + p\frac{\partial H_c}{\partial s} = - kx - \gamma \, p}\\ {\dot{s}} & =& H_c-p {\displaystyle \frac{\partial H_c}{\partial p} = \frac{m_{0}^{2}c^2}{\sqrt{p^2 +m_0^2c^2}} - \frac{1}{2}kx^2} \end{array}\right. \qquad . \end{aligned}$$

(63)

Using an elementary algebra we obtain

$$\begin{aligned} 1 - \frac{{\dot{x}}^2}{c^2} = \frac{m_{0}^{2}c^2}{p^2 + m_{0}^{2}c^2}. \end{aligned}$$

Then taking into account the first two equations of the set (63) we obtain

$$\begin{aligned} \ddot{x}=\left( \frac{c}{\sqrt{p^2 +m_0^2c^2}}- \frac{cp^2}{(p^2 +m_0^2c^2)^{3/2}}\right) \dot{p}=\frac{m_0^2c^3}{(p^2 +m_0^2c^2)^{3/2}}(-kx-\gamma p), \end{aligned}$$

which can be expressed as

$$\begin{aligned} \ddot{x} + \gamma {\dot{x}}\left( 1 - \frac{{\dot{x}}^2}{c^2} \right) + \left( 1 - \frac{{\dot{x}}^2}{c^2} \right) ^{3/2}\omega ^2 x = 0, \quad \omega ^2 = \frac{k}{m_0}. \end{aligned}$$

(64)

We recover the non-relativistic limit from the limit $\dfrac{{\dot{x}}}{c} \rightarrow 0$, but actually we have being working here in the semi-relativistic framework, i.e., by keeping the relativistic corrections up to order $1/c^2$. We need general relativity to properly formulate the relativistic harmonic oscillator, where we will deal curved space-time metric rather than flat metric (Chanda and Guha 2018).

3 Levinson–Smith Equation and Contact Hamiltonian Formalism

As indicated above, the Levinson–Smith Eq. (3) is a generalization of the Liénard Eq. (1) by containing a term proportional to ${\dot{x}}^2$. The damped position-dependent mass Eq. (51) is a typical example of the Levinson–Smith equation. In general this equation cannot be framed in terms of the standard (contact) Hamiltonian formulation, but it is rather connected to a non-standard Hamiltonian formulation. As mentioned in Mitra et al. (2024) the particular Levinson–Smith Eq. (4) given by

$$\begin{aligned} \ddot{x} + p(x){\dot{x}}^{2} + q(x){\dot{x}} + g(x) = 0, \end{aligned}$$

(65)

may be expressed as the planar system

$$\begin{aligned} \left\{ \begin{array}{rcl} {\dot{x}}& =& y \\ {\dot{y}} & =& - p(x)y^2 - q(x) y - g(x) \end{array}\right. \qquad , \end{aligned}$$

(66)

with an associated vector field in ${\mathbb {R}}^2$

$$\begin{aligned} \Gamma _{\textrm{LS}} = y\frac{\partial }{\partial x} - (p(x)y^2 + q(x)y + g(x))\frac{\partial }{\partial y}. \end{aligned}$$

(67)

It may be convenient to replace the system (66) by an extended one involving the time as a new coordinate

$$\begin{aligned} \left\{ \begin{array}{rcl} \dfrac{dt}{d\tau }& =& 1\\ \dfrac{dx}{d\tau }& =& y \\ \dfrac{dy}{d\tau }& =& - p(x)y^2 - q(x) y - g(x) \end{array}\right. \end{aligned}$$

(68)

with an associated vector field in ${\mathbb {R}}^3$

$$\begin{aligned} {{\bar{\Gamma }}}_{\textrm{LS}} =\frac{\partial }{\partial t}+ y\frac{\partial }{\partial x} - (p(x)y^2 + q(x)y + g(x))\frac{\partial }{\partial y}. \end{aligned}$$

(69)

If we implicitly assume the volume form in ${\mathbb {R}}^3$, $\omega =dt\wedge dx\wedge dy$, the planar system given by (66) admits a Jacobi Last Multiplier (JLM), ${\mathcal {M}}$, which, by definition, is a solution to the equation (Jacobi 1845, 2009)

$$\begin{aligned} {{\bar{\Gamma }}}_{\textrm{LS}}(\log {\mathcal {M}})+\mathrm{div\,} {{\bar{\Gamma }}}_{\textrm{LS}}=0. \end{aligned}$$

(70)

Note that this equation shows that another function ${\mathcal {M}} '$ is also a JLM if and only if the quotient ${\mathcal {M}}'/ {\mathcal {M}}$ is a first integral for ${{\bar{\Gamma }}}_{\textrm{LS}}$.

For instance, for the classical damped harmonic oscillator, $\ddot{x}+\lambda \, \dot{x}+\omega ^2\, x=0$, with $\lambda \in {{\mathbb {R}}}, 0<\omega \in {\mathbb {R}}$, corresponding $q(x)=\lambda , \, p(x)=0, \, g(x)=\omega ^2\, x$, it is described by the vector field

$$\begin{aligned} {{\bar{\Gamma }}}_{\textrm{do}} =\frac{\partial }{\partial t}+y\frac{\partial }{\partial x}-(\lambda \, y+\omega ^2\, x)\frac{\partial }{\partial y}, \end{aligned}$$

(71)

hence we have that $\mathrm{div\,} {{\bar{\Gamma }}}_{\textrm{do}} =-\lambda$ and then in this case there is a JLM given by ${\mathcal {M}}=e^{\lambda t}$.

The more general case of a mechanical type system with an additional damping term, i.e., the second-order differential equation

$$\begin{aligned} \ddot{x}+\lambda \, \dot{x}+p(x)\,\dot{x}^2+g(x)=0, \end{aligned}$$

(72)

where p is a positive function, can be dealt with in a similar way. This is a particular case of (65) with a function $q(x)=\lambda \in {\mathbb {R}}$. Now the relevant vector field is

$$\begin{aligned} {{\bar{\Gamma }}}_{\textrm{dmt}} =\frac{\partial }{\partial t}+y\frac{\partial }{\partial x}-(\lambda \, y+p(x)y^2+g(x))\frac{\partial }{\partial y}, \end{aligned}$$

(73)

for which $\mathrm{div\,} {{\bar{\Gamma }}}_{\textrm{dmt}} =-2p(x)y-\lambda$. In the particular case $\lambda =0$, ${{\bar{\Gamma }}}_{\textrm{dmt}}$ reduces to

$$\begin{aligned} {{\bar{\Gamma }}}_{\textrm{dmt}}^0=\frac{\partial }{\partial t}+y\frac{\partial }{\partial x}-\left( p(x)\, y^2+g(x)\right) \frac{\partial }{\partial y}, \end{aligned}$$

and then, if ${\mathbb {R}}^3$ is endowed with the volume form $\omega =dt\wedge dx\wedge dy$, as $\mathrm{div\,}\Gamma _{\textrm{dmt}}^0=- 2p(x)\, y$, the function $R_0$ of the variable x defined by

$$\begin{aligned} R_0(x)=\exp \left( 2\int ^xp(\zeta )\, d\zeta \right) , \end{aligned}$$

is such that

$$\begin{aligned} {{\bar{\Gamma }}}_{\textrm{dmt}}^0(R_0)= 2p(x)\exp \left( 2\int ^xp(\zeta )\, d\zeta \right) \, y, \end{aligned}$$

and as ${{\bar{\Gamma }}}_{\textrm{dmt}}^0 (R_0)+R_0\, \mathrm{div\,}\Gamma _{\textrm{dmt}}^0=0$, the positive function $R_0$ is a Jacobi multiplier for the vector field ${{\bar{\Gamma }}}_{\textrm{dmt}}^0$.

Coming back to the complete Eq. (72) we will have a vector field ${{\bar{\Gamma }}}_{\textrm{dmt}}\in {{\mathfrak {X}}}({\mathbb {R}}^3)$ given by (73). We can look for a JLM for ${{\bar{\Gamma }}}_{\textrm{dmt}}$ of the form $R(x,t)=u(t) R_0(x)$ and then, as

$$\begin{aligned} \mathrm{div\,}{{\bar{\Gamma }}}_{\textrm{dmt}}(t,x,v)=-\lambda - 2p(x)\, y,\qquad {{\bar{\Gamma }}}_{\textrm{dmt}} (R)=\dot{u}\, R_0+u\, y\,\frac{dR_0}{dx}=\dot{u}\, R_0+2u\, y\,p\,R_0, \end{aligned}$$

from where we see that $R\,\mathrm{div\,}{{\bar{\Gamma }}}_{\textrm{dmt}} +{{\bar{\Gamma }}}_{\textrm{dmt}}(R)=(-\lambda u+\dot{u})R_0$, and therefore if the function u(t) is of the form $u(t)=e^{\lambda t}$, then the function $R=e^{\lambda t}R_0$ is a Jacobi multiplier for the vector field ${{\bar{\Gamma }}}_{\textrm{dmt}}$ with respect to the volume form $dt\wedge dx \wedge dv$.

The particular case of (65) with $p(x)\equiv 0$, has been studied in Cariñena and Guha (2019) and it has been shown that if the Chiellini condition holds for a real number $\alpha$:

$$\begin{aligned} \frac{d}{dx}\left( \frac{g(x)}{q(x)}\right) =-\alpha \,(\alpha +1)q(x)\,\quad -1\ne \alpha \ne 0, \end{aligned}$$

(74)

then there is a JLM of the form ${\mathcal {M}}=(y-W(x))^{1/\alpha }$, where $W(x)=\dfrac{1}{\alpha }\,\dfrac{g(x)}{q(x)}$. Note also that if $-\alpha (\alpha +1)=r$, when. $r<1/4$ exists another real number $\beta$ such that $-\beta (\beta +1)=r$ and therefore another similar JLM, therefore there is a first integral given by the quotient of both JLMs.

Note also that, as indicated in Cariñena et al. (2023), an Eq. (65) with $p(x)\ne 0$, recently studied in Mitra et al. (2024), can be related to a similar one with $p(x)\equiv 0$ by means of a pure Sundman transformation, a kind of infinitesimal time-reparametrization (see e.g., Cariñena et al. (2023, 2022)). Actually, if we define a Sundman transformation by

$$\begin{aligned} \frac{d\tau }{dt}=\exp \left( - \int ^x p(\zeta )\, d\zeta \right) =e^{-P(x)}, \end{aligned}$$

(75)

then,

$$\begin{aligned} \frac{dx}{dt}=e^{-P(x)} \frac{dx}{d\tau }, \qquad \frac{d^2x}{dt^2}=e^{-P(x)}\frac{d}{d\tau }\left( e^{-P(x)}\frac{dx}{d\tau }\right) =e^{-2P(x)}\left( \frac{d^2x}{d\tau ^2}-p\,\left( \frac{dx}{d\tau }\right) ^2\right) , \end{aligned}$$

and consequently the new equation is

$$\begin{aligned} \frac{d^2x}{dt^2} + p(x)\left( \frac{dx}{dt}\right) ^2 + q(x)\frac{dx}{dt}+ g(x) =e^{-2P(x)}\frac{d^2x}{d\tau ^2}+e^{-P(x)} q(x)\, \frac{dx}{d\tau }+g(x)=0, \end{aligned}$$

which looks as the usual Liénard equation. The Chiellini condition for this equation is:

$$\begin{aligned} \frac{d}{dx}\left( e^{P(x)}\frac{g(x)}{q(x)}\right) =-\alpha (\alpha +1)e^{P(x)}q(x), \end{aligned}$$

(76)

and if this condition (76) is satisfied, the JLM will be ${\mathcal {M}}=(e^{P(x)}y-W(x))^{1/\alpha }$, where

$$\begin{aligned} W(x)=\dfrac{1}{\alpha }e^{P(x)} \dfrac{g(x)}{q(x) }, \end{aligned}$$

and the JLM is then

$$\begin{aligned} {\mathcal {M}}=e^{(1/\alpha )P(x)}\left( y-\frac{1}{\alpha }\frac{g(x)}{q(x)}\right) ^{1/\alpha }. \end{aligned}$$

Once again, as above, if $-\alpha (\alpha +1)=r$, when $r<1/4$, there exists another real number $\beta$ such that $-\beta (\beta +1)=r$ and then another similar JLM, therefore there is a first integral given by the quotient of both JLMs (Mitra et al. 2024).

Proposition 4

The Levinson–Smith equation

$$\begin{aligned} \ddot{x} + p(x){\dot{x}}^{2} + q(x){\dot{x}} + g(x) = 0, \end{aligned}$$

satisfying Chiellini integrability condition (76) admits a Lagrangian description by a non-standard Lagrangian function of the form

$$\begin{aligned} L(x,\dot{x}) = \frac{1}{((1/\alpha )+1)((1/\alpha )+2)} e^{(1/\alpha ) P(x)} \left( {\dot{x}} - \frac{1}{\alpha }\frac{g(x)}{q(x)}\right) ^{(1/\alpha )+2}, \end{aligned}$$

(77)

where $P(x) = {\displaystyle \int ^x} p(x^{\prime })dx^{\prime }$.

Proof

The proof is based on the fact that for any second-order ordinary differential equation of the form $\ddot{x}={\mathcal {F}}(x, {\dot{x}})$ which admits a JLM, ${\mathcal {M}}$, the latter is connected to a Lagrangian for the system, uniquely defined up to addition of a gauge term, via the relation ${\mathcal {M}}(x,y)=\partial ^2\,L/\partial y^2$ (see (Nucci and Leach 2008, 2009; Nucci and Tamizhmani 2010) or Cariñena and Fernández-Núñez (2021) for a more recent version). Consequently (see e.g., Cariñena and Fernández-Núñez (2021) and references therein) a Lagrangian L is given, up to addition of a gauge term, by

$$\begin{aligned} L(x,y)=\int ^y d\zeta \int ^\zeta {\mathcal {M}}(x,\eta )\, d\eta -V(x), \end{aligned}$$

and V(x) is uniquely defined, up to addition of a constant. In the present case it follows therefore that

$$\begin{aligned} L(x,y)=\frac{1}{((1/\alpha )+1)((1/\alpha )+2)} e^{(1/\alpha ) P(x)} \left( y - \frac{1}{\alpha }\frac{g(x)}{q(x)}\right) ^{(1/\alpha )+2}-V(x). \end{aligned}$$

To determine the function V(x) we substitute this form of the Lagrangian into the Euler-Lagrange equation and by comparing the result with the given equation, i.e., (65), the value of V(x) can be determined. In our case it turns out that $V(x)=0$ and hence the proof. $\square$

This is not a natural Lagrangian or, in other words, this is a non-standard Lagrangian. The corresponding conserved energy function $E_L = y \dfrac{\partial L}{\partial y} - L$ is given by

$$\begin{aligned} E_L(x,y) = \frac{e^{(1/\alpha ) P(x)}}{(1/\alpha ) + 2} \left( y - \frac{1}{\alpha }\frac{g(x)}{q(x)}\right) ^{(1/\alpha )+1} \Big ( y + \frac{1}{\alpha ((1/\alpha ) + 1)}\frac{g(x)}{q(x)}\Big ), \end{aligned}$$

(78)

which yields

$$\begin{aligned} E_L & = \frac{e^{(1/\alpha ) P(x)}}{(\alpha +1)(1/\alpha ) + 2} \left( y -\frac{1}{\alpha }\frac{g(x)}{q(x)}\right) ^{1/\alpha }\nonumber \\ & \quad \ s\left( y^2 - \frac{1}{\alpha (\alpha + 1)}\frac{g(x)}{q(x)}y - \frac{1}{\alpha (\alpha + 1)} \left( \frac{g(x)}{q(x)}\right) ^2 \right) . \end{aligned}$$

(79)

Such energy function $E_L$ is an integral of motion of the Levinson–Smith equation.

We may recover the energy function of the Liénard equation (see (17) of Cariñena and Guha (2019)) from the reduction of the energy function of the Levinson–Smith equation.

Corollary 1

When $p(x) = 0$, the energy $E_L$ reduces to the energy function of the Liénard equation

$$\begin{aligned} E_{\textrm{Lienard}} = \frac{1}{(1/\alpha ) + 2} \left( y - \frac{1}{\alpha }\frac{g(x)}{q(x)}\right) ^{(1/\alpha )+1} \Big ( y + \frac{1}{\alpha ((1/\alpha ) + 1)}\frac{g(x)}{q(x)}\Big ). \end{aligned}$$

(80)

The expression for the Hamiltonian can now be easily derived via the Legendre transformation, i.e.,

$$\begin{aligned} \pi = \frac{e^{(1/\alpha ) P(x)}}{(1/\alpha ) + 1} \left( y - \frac{1}{\alpha }\frac{g(x)}{q(x)}\right) ^{(1/\alpha )+1}, \end{aligned}$$

which can be inverted as

$$\begin{aligned} y=\left[ \left( (1/\alpha ) + 1\right) e^{-(1/\alpha ) P(x)}\pi \right] ^{\frac{\alpha }{\alpha +1}}+\frac{1}{\alpha }\frac{g(x)}{q(x)}, \end{aligned}$$

and then the energy has the following appearance

$$\begin{aligned} H=\frac{\alpha +1}{2\alpha +1}\pi \left( y- \frac{1}{\alpha }\frac{g(x)}{q(x)}+\frac{2\alpha +1}{\alpha (\alpha +1)} \frac{g(x)}{q(x)}\right) , \end{aligned}$$

and hence

$$\begin{aligned} H=\frac{\alpha +1}{2\alpha +1}\pi \left( \left[ \left( (1/\alpha ) + 2\right) e^{-(1/\alpha ) P(x)}\pi \right] ^{\frac{\alpha }{\alpha +1}} +\frac{2\alpha +1}{\alpha (\alpha +1)} \frac{g(x)}{q(x)}\right) . \end{aligned}$$

3.1 Contact Hamiltonians Formulation of the Levinson–Smith Equation and Reductions

Let us consider the Levinson–Smith equation with $p(x) = k$ (say) satisfying Chiellini condition. The quadratic part of its energy function $E_L$ yields the following function

$$\begin{aligned} H(x,y)= e^{kx/\alpha }\left( y^2 - \frac{1}{\alpha (\alpha + 1)}y\frac{g(x)}{q(x)} - \frac{1}{\alpha (\alpha + 1)}\left( \frac{g(x)}{q(x)}\right) ^2 \right) , \end{aligned}$$

(81)

where we ignore the other irrelevant factors, these are constant terms which do not contribute in the dynamics of motion.

We can introduce a new variable ${\bar{s}}$, called entropy, and consider in $\mathbb {R}^3$ the surface defined by

$$\begin{aligned} {\bar{s}}(x,y) = ye^{kx/\alpha }\frac{g(x)}{q(x)}. \end{aligned}$$

(82)

Thus the contact function corresponding to (81) is given as

$$\begin{aligned} H_{\textrm{CLS}} (x,y,{\bar{s}})= e^{kx/\alpha }\left( y^2 + \gamma ^2 \left( \frac{g(x)}{q(x)}\right) ^2 \right) + \gamma ^2\,{\bar{s}}, \end{aligned}$$

(83)

where $\gamma$ and $\alpha$ are related by $(\gamma ^2)^{-1} = - \alpha (\alpha + 1)$. This function is such that the pull-back of the function $H_{\textrm{CLS}} (x,y,{{\bar{s}}})$ on the surface defined by (82) coincides with H(x, y).

3.1.1 Illustration: Generalized Damped Oscillator

Case 1.  Consider the particular case $g(x) = x$ and $q(x) = \gamma$, and define $k'=k/\alpha$, hence the Hamiltonian (81) becomes

$$\begin{aligned} H_{1} = e^{k'x}\left( y^2 + \gamma \,xy +x^2 \right) , \end{aligned}$$

(84)

using the entropy function $s:=\gamma \, {\bar{s}}= xy \,e^{k'x}$, which yields the conformal Hamiltonian partner

$$\begin{aligned} H_{\textrm{CLS}} = e^{k'x}\left( y^2 + x^2 \right) + \gamma \, s. \end{aligned}$$

(85)

The system of differential equations for the integral curves of the contact Hamiltonian system with the contact 1-form $\eta =ds+y\, dx$, is given by

$$\begin{aligned} \left\{ \begin{array}{rcl} {\dot{x}} & =& 2e^{k'x}y\\ {\dot{y}} & =& - k' e^{k'x}\big (y^2 + x^2) - 2xe^{k'x} - \gamma \, y\\ {\dot{s}} & =& e^{k'x}( y^2 - x^2 - \gamma \,s). \end{array}\right. \quad . \end{aligned}$$

(86)

Using the first two equations we obtain

$$\begin{aligned} \ddot{x} - \frac{1}{2}{\dot{x}}^{2} + \gamma \, {\dot{x}} + 2k'e^{2k'x}x^2 + 4xe^{2k'x} = 0. \end{aligned}$$

(87)

This reduces to the standard damped oscillator equation for $k = 0$ with two additional scalings, namely, $\gamma \mapsto 2\gamma$ and $t \mapsto t/2$.

Case 2.  Let us take $p(x) = \dfrac{m^{\prime }(x)}{m(x)}$. This immediately yields

$$\begin{aligned} P(x)=\exp \left( \int ^x\frac{m^{\prime }(\zeta )}{m(\zeta )}d\zeta \right) =m(x), \end{aligned}$$

and the integration constant can be chosen to cancel the factor $\alpha$, i.e., we can replace $P(x)/\alpha$ by m(x). Thus the contact Hamiltonian obtained from $H = m(x)(y^2 + \gamma xy + x^2)$ is given by

$$\begin{aligned} H_{\textrm{CPM}} = m(x) \big ( y^2 + x^2 \big ) + \gamma \,s, \quad s = m(x)xy, \end{aligned}$$

(88)

which yields the following system of equations of motion

$$\begin{aligned} \left\{ \begin{array}{rcl} {\dot{x}} & =& 2m(x)y\\ {\dot{y}} & =& - m^{\prime }(x) (y^2 + x^2) - 2m(x)x - \gamma \, y\\ {\dot{s}} & =& m(x)\big ( y^2 - x^2) - \gamma \, s. \end{array}\right. \end{aligned}$$

(89)

This can be reduced to a position-dependent mass equation

$$\begin{aligned} \ddot{x} - \frac{1}{2}\frac{m^{\prime }(x)}{m(x)}{\dot{x}}^{2} + \gamma {\dot{x}} + 2m(x)m^{\prime }(x)x^2 + 4m^2(x)x = 0. \end{aligned}$$

(90)

4 GENERIC Formalism of Levinson–Smith Equation

It will be shown in this section that this form of Liénard equation can be realized through a mixture of a Hamiltonian and a gradient flow. Our discussion is based on the GENERIC (General Equation for the Non-Equilibrium Reversible-Irreversible Coupling) formulation of time evolution equations manufactured for non-equilibrium systems (Öttinger and Grmela 1997a, b; Öttinger 2018; Mielke 2011; Mielke et al. 2024). As the geometric approach to this formulation is not well-known we first recall the main ingredients.

In the geometric approach to dynamics in classical mechanics the Hamiltonian formulation in symplectic manifolds $(M,\omega )$ is the most frequently used but its generalization to the framework of Poisson manifolds $( M,\Lambda )$ is also very useful. The symplectic 2.form, i.e., a non-degenerate closed 2-form $\omega$, is replaced by a skew-symmetric 2-times contravariant tensor field, $\Lambda$, and the condition replacing closedness of $\omega$ is the vanishing of the Schouten bracket, $[\Lambda ,\Lambda ]_{\textrm{Sch}}=0$, for the bivector field $\Lambda$. Nondegeneracy is lost in this Poisson approach, but it is also possible to introduce a Poisson bracket by $\{F,G\}=\Lambda (dF,dG)$ which may admit Casimir functions ${\mathcal {C}}$ such that $\{{\mathcal {C}},G\}=0, \forall G\in C^\infty (M)$. Jacobi identity holds as a consequence of the $[\Lambda ,\Lambda ]_{\textrm{Sch}}=0$ condition. The dynamical vector field defined by the Hamiltonian function $H\in C^\infty (M)$ is given by $X_H(G)= \Lambda (dG,dH),\forall G\in C^\infty (M)$. This system is conservative in the sense that $X_H(H) = \{H,H\}=0$, by skew-symmetry of the Poisson bracket.

If the local expression of the Poisson bivector field $\Lambda$ is

$$\begin{aligned} \Lambda =L^{ij}(x)\frac{\partial }{\partial x^i}\wedge \frac{\partial }{\partial x^j},\quad L^{ij} +L^{ji}=0, \end{aligned}$$

then that of the vanishing of the Schouten bracket, $[\Lambda ,\Lambda ]_{\textrm{Sch}}=0$, is

$$\begin{aligned} [\Lambda ,\Lambda ]_{\textrm{Sch}}=L^{ik}{\displaystyle \frac{\partial L^{lj}}{\partial x^k}+L^{jk}\frac{\partial L^{il}}{\partial x^k}+L^{kl}\frac{\partial L^{ij}}{\partial x^k}}=0, \end{aligned}$$

that of the Poisson bracket is

$$\begin{aligned} \{F,G\}=L^{ij}(x)\frac{\partial F}{\partial x^i}\, \frac{\partial G}{\partial x^j}, \end{aligned}$$

and that of the Hamiltonian vector field $X_H$ is

$$\begin{aligned} X_H=L^{ij}(x)\frac{\partial H}{\partial x^j}\, \frac{\partial }{\partial x^i}. \end{aligned}$$

(91)

Its integral curves are the solutions of

$$\begin{aligned} \frac{dx^i}{dt}=L^{ij}(x)\frac{\partial H}{\partial x^j},\quad i=1,\ldots ,n, \end{aligned}$$

(92)

and the evolution of the function $G \in C^\infty (M)$ along such integral curves is given by

$$\begin{aligned} \frac{dG}{dt}=L^{ij}(x)\frac{\partial G}{\partial x^i}\, \frac{\partial H}{\partial x^j}=\{G,H\}. \end{aligned}$$

(93)

In order to have room for the evolution of nonconservative systems, as it is the case of thermodynamics, we must replace the Hamiltonian vector field by a different one, for instance adding a new dissipating vector field. This is what it is made in the GENERIC (general equation for the non-equilibrium reversible-irreversible coupling) formulation of the dynamics (Esen et al. 2022a, b; Öttinger and Grmela 1997a, b; Öttinger 2018). It was shown in Kaufman (1984); Morrison (2009, 1986), Guha (2007) that the appropriate geometric framework for studying dissipation is that of pseudo-Riemannian manifolds and the use of a dissipative bracket (Kaufman 1984).

Recall that in a Riemannian manifold (M, g) the so-called gradient vector fields are very relevant. They are defined by the relation

$$\begin{aligned} g(\mathrm{grad\,} F,\cdot )=dF(\cdot ). \end{aligned}$$

(94)

We can forget the nondegeneracy condition for g and even go to the dual formulation and define a symmetric two-times contravariant positive tensor field, ${\mathcal {G}}$, and define the dissipative vector field $Y_F$ defined by the function $F\in C^\infty (M)$ (they are called metric systems in Morrison (1986)) by means of

$$\begin{aligned} Y_F(G)= {\mathcal {G}}(dF,dG),\quad \forall G\in C^\infty (M). \end{aligned}$$

(95)

If the local expression of the tensor field ${\mathcal {G}}$ is

$$\begin{aligned} {\mathcal {G}}=\textrm{M}^{ij} (x)\,\frac{\partial }{\partial x^i}\otimes \frac{\partial }{\partial x^j},\quad \textrm{M}^{ij} =\textrm{M} ^{ji}, \end{aligned}$$

(96)

then the local expression of the vector field $Y_F$ is

$$\begin{aligned} Y_F= \textrm{M}^{ij}(x)\frac{\partial F}{\partial x^i}\,\frac{\partial }{\partial x^j}, \end{aligned}$$

(97)

and it is possible to define a symmetric pairing in $C^\infty (M)$, called dissipative bracket, $\langle \langle F,G\rangle \rangle = {\mathcal {G}}(dF,dG)$, with a coordinate expression

$$\begin{aligned} \langle \langle F,G\rangle \rangle =\textrm{M}^{ij}(x)\frac{\partial F}{\partial x^i}\, \frac{\partial G}{\partial x^j}= \langle \langle G,F\rangle \rangle , \end{aligned}$$

(98)

that, because of the positivity assumption on ${\mathcal {G}}$, implies that

$$\begin{aligned} \langle \langle F,F\rangle \rangle \ge 0,\quad \forall F \in C^\infty (M). \end{aligned}$$

(99)

In GENERIC formulation it is assumed that two functions E and S (recalling the rôle of the energy and the entropy) are given and then the dynamical vector field is $X_E+Y_S=\{\cdot ,E\}+ \langle \langle \cdot ,S\rangle \rangle$, whose integral curves are solutions of

$$\begin{aligned} \frac{dx^i}{dt}=L^{ij}(x)\frac{\partial E}{\partial x^j}+\textrm{M}^{ij}(x)\frac{\partial S}{\partial x^j}. \end{aligned}$$

(100)

It was also originally postulated, on physical grounds, that both functions E and S are such that $\{E,S\}=0$, and later on this condition was extended to $\{F,S\}=0, \forall F\in C^\infty (M)$, or, in other words, the assumption is that the entropy function S is a Casimir function, namely $\Lambda (dS,\cdot )=0$, or in coordinates,

$$\begin{aligned} L^{ij}(x)\frac{\partial S}{\partial x^j}=0. \end{aligned}$$

(101)

Finally, in the GENERIC formulation the additional constraint on the function E is introduced: $\langle \langle E,F\rangle \rangle =0, \,\forall F\in C^\infty (M)$. In local coordinates,

$$\begin{aligned} \textrm{M}^{ij}(x)\frac{\partial E}{\partial x^i}=0. \end{aligned}$$

(102)

Equations (101) and (102) indicate the conservation of the entropy by the reversible dynamics and the conservation of the total energy in a closed system by the irreversible dynamics, respectively.

The particular example of a one-dimensional harmonic oscillator in presence of friction,

$$\begin{aligned} \left\{ \begin{array}{rcl}\dfrac{dq}{dt}& =& p\\ \dfrac{dp}{dt}& =& -q-\gamma \, p\\ \dfrac{ds}{dt}& =& \gamma \, p^2\end{array} \right. \end{aligned}$$

(103)

was presented in Öttinger (2018) with $\Lambda ={\displaystyle \frac{\partial }{\partial q}\wedge \frac{\partial }{\partial p}}$,

$$\begin{aligned} {\mathcal {G}}= \gamma \left( \frac{\partial }{\partial p}\otimes \frac{\partial }{\partial p}-\left( \frac{\partial }{\partial p}\otimes \frac{\partial }{\partial s}+ \frac{\partial }{\partial s}\otimes \frac{\partial }{\partial p}\right) + p^2\frac{\partial }{\partial s}\otimes \frac{\partial }{\partial s}\right) , \end{aligned}$$

and

$$\begin{aligned} E(q,p,s)=\frac{p^2}{2} +\frac{1}{2} q^2+s, \qquad S(q,p,s)=s, \end{aligned}$$

which satisfy (101) and (102), and we obtain the vector fields

$$\begin{aligned} X_E= p\frac{\partial }{\partial q}-q\frac{\partial }{\partial p},\qquad Y_S=-\gamma p\frac{\partial }{\partial p}+\gamma \, p^2\frac{\partial }{\partial s}, \end{aligned}$$

and then one can check that actually the integral curves of the vector field $X_F+Y_S$ are the solutions of (103).

4.1 Illustration: Construction of GENERIC Formulation for Liénard Equation

The similarity of Liénard equation with damped harmonic oscillator equation suggests us to search for a GENERIC formulation for Liénard equation.

Given a function $f\in C^\infty (Q)$ we can consider the Poisson bivector field on $TQ\times {\mathbb {R}}$ that in a coordinate system induced from one on Q is $\Lambda ={\displaystyle \frac{\partial }{\partial x}\wedge \frac{\partial }{\partial y}}$ with

$$\begin{aligned} {\mathcal {G}}= f(x)\left( \frac{\partial }{\partial y}\otimes \frac{\partial }{\partial y}-m(x)y\left( \frac{\partial }{\partial y}\otimes \frac{\partial }{\partial s}+ \frac{\partial }{\partial s}\otimes \frac{\partial }{\partial y}\right) + y^2\frac{\partial }{\partial s}\otimes \frac{\partial }{\partial s}\right) ,\nonumber \\ \end{aligned}$$

(104)

and

$$\begin{aligned} E(x,y,s)=\frac{1}{2} \left( m(x) y^2-\frac{1}{\alpha (\alpha +1)}\left( \frac{g}{f}\right) ^2 \right) +s, \qquad S(x,y,s)=s, \end{aligned}$$

(105)

where it is assumed the Chiellini condition

$$\begin{aligned} \frac{d}{dx}\left( \frac{g}{f}\right) =-\alpha (\alpha +1) f. \end{aligned}$$

(106)

One can readily check that the conditions (101) and (102) are satisfied, and we obtain the following coordinate expression for the vector fields $X_E$ and $Y_S$:

$$\begin{aligned} & X_E= m(x)y \frac{\partial }{\partial x}-\left( \frac{1}{2} m'(x)y+g(x)\right) \frac{\partial }{\partial y},\\ & Y_S=-m(x) f(x)y\frac{\partial }{\partial y}+f(x)\,m^2(x) y^2\frac{\partial }{\partial s}, \end{aligned}$$

and therefore the integral curves of the dynamical vector field $X_E+Y_S$ are solutions of the system

$$\begin{aligned} \left\{ \begin{array}{rcl} \dfrac{dx}{dt}& =& m(x)\, y\\ \dfrac{dy}{dt}& =& - \dfrac{1}{2} m'(x)y^2-g(x) \\ \dfrac{ds}{dt}& =& f(x)\,m^2(x) y^2 \end{array} \right. \quad . \end{aligned}$$

(107)

So, we have obtained the following result:

Proposition 5

The GENERIC equations of motion defined by the Poisson bivector field $\Lambda$ on $TQ\times {\mathbb {R}}$ and the symmetric two-times contravariant tensor field ${\mathcal {G}}$ defined in a coordinate system in $TQ\times {\mathbb {R}}$ induced from one on Q by $\Lambda ={\displaystyle \frac{\partial }{\partial x}\wedge \frac{\partial }{\partial y}}$ and (104), together with the energy and entropy functions (105) yield, when Chiellini condition (106) holds, to the system (107) of a Levinson–Smith equation

Note that from the two first equations of system (107) we see that the projections of solutions of such system on the base manifold Q are the solutions of the second-order differential equation

$$\begin{aligned} \ddot{x} + \frac{1}{2}\frac{m^{\prime }(x)}{m(x)}{\dot{x}}^2 + f(x){\dot{x}} + g(x) = 0, \end{aligned}$$

(108)

where the coefficient of the ${\dot{x}}^{2}$ term is $p(x) = \dfrac{1}{2}\dfrac{m^{\prime }(x)}{m(x)}$.

We recover the damped oscillator GENERIC construction of the damped oscillator $\ddot{x} + \gamma \, {\dot{x}} + x = 0$ given in Öttinger (2018) for $f(x) = \gamma$, $m(x) =1$ and $g(x) = x$.

Example 1

A position-dependent mass system and the Levinson–Smith Equation.

A one-dimensional position-dependent mass system is described in the Hamiltonian formulation by a natural Hamiltonian that in coordinates of a chart on $T^*Q$ induced from a chart on Q is

$$\begin{aligned} H(x, p)=\frac{p^2}{2m(x)}+V(x), \end{aligned}$$

and if we consider the Poisson bivector field $\Lambda ={\displaystyle \frac{\partial }{\partial q}\wedge \frac{\partial }{\partial p}}$ on $T^*Q$ we will have the associated Hamilton equations

$$\begin{aligned} \left\{ \begin{array}{rcl} \dfrac{dx}{dt}& =& \dfrac{p}{m(x)}\\ \dfrac{dp}{dt}& =& \dfrac{m^\prime (x)}{2m^2(x)}p^2-V^\prime (x) \end{array}\right. \quad , \end{aligned}$$

(109)

and then the projections on the base manifold Q of the curves solution of such system satisfy the second-order differential equation

$$\begin{aligned} \ddot{x}+\frac{1}{2} \frac{m'(x)}{m(x)}\dot{x}^2+\frac{1}{m(x)}V'(x)=0. \end{aligned}$$

Of course, if the potential function V is such that $V'/m=g$, i.e., $V(x)={\displaystyle \int ^x m(z)g(z)\, dz}$, such equation is a particular case of (4).

It is also possible to study the generic formulation of the case of an additional dissipative term in the equation,

$$\begin{aligned} \ddot{x}+\frac{1}{2} \frac{m'(x)}{m(x)}\dot{x}^2+f(x)\dot{x}+\frac{1}{m(x)}V'(x)=0. \end{aligned}$$

(110)

We can still choose the (extension to $T^*Q\times {{\mathbb {R}}}$ of the) Poisson bivector field $\Lambda ={\displaystyle \frac{\partial }{\partial q}\wedge \frac{\partial }{\partial p}}$ and functions E and S on $T^*Q\times {{\mathbb {R}}}$ of the form

$$\begin{aligned} E(x,p,s)=\frac{p^2}{2m(x)}+V(x)+U(s),\qquad S(x,p,s)=s, \end{aligned}$$

(111)

and we can check that there exists a choice for the function U(s) such that the integral curves of the dynamical vector field $X_E+Y_S$ corresponding to

$$\begin{aligned} {\mathcal {G}} & = f(x)m(x)\left( (U'(s))^2 \frac{\partial }{\partial p}\otimes \frac{\partial }{\partial p}-\frac{pU'(s)}{m(x)}\left( \frac{\partial }{\partial p}\otimes \frac{\partial }{\partial s}+ \frac{\partial }{\partial s}\otimes \frac{\partial }{\partial p}\right) \right. \nonumber \\ & \left. +\frac{p^2}{m^2(x)} \frac{\partial }{\partial s}\otimes \frac{\partial }{\partial s}\right) , \end{aligned}$$

(112)

are solutions of the second-order differential equation for the damped position-dependent vector field (110), because then

$$\begin{aligned} & X_E=\frac{p}{m(x)}\frac{\partial }{\partial x}+ \left( \frac{m^\prime (x)}{2m^2(x)}p^2-V^\prime (x)\right) \frac{\partial }{\partial p},\\ & Y_S= f(x)m(x)\left( -\frac{pU'(s)}{m(x)}\frac{\partial }{\partial p}+\frac{p^2}{m^2(x)}\frac{\partial }{\partial s}\right) , \end{aligned}$$

and the integral curves of the vector field $X_E+Y_S$ are solutions of

$$\begin{aligned} \left\{ \begin{array}{rcl} \dfrac{dx}{dt}& =& \dfrac{p}{m(x)}\\ \dfrac{dp}{dt}& =& \dfrac{m^\prime (x)}{2m^2(x)}p^2-V^\prime (x) -f(x)\,p\,U'(s)\\ \dfrac{ds}{dt}& =& f(x)\dfrac{p^2}{m(x)} \end{array}\right. \quad , \end{aligned}$$

(113)

and then the projections on the base manifold Q of the curves solution of such system satisfy the second-order differential equation

$$\begin{aligned} \ddot{x}+\frac{1}{2} \frac{m'(x)}{m(x)}\dot{x}^2+\frac{1}{m(x)}V'(x)+f(x)\, U'(s)\,\dot{x}=0, \end{aligned}$$

which shows that the choice $U'(s)=1$ leads to (110).

5 Discussion

This is a continuation of our program started in Cariñena and Guha (2019, 2024), where the geometry of the Liénard equation was studied. In this article we studied the Levinson–Smith equation. This equation is related to the Liénard-type second-order nonlinear differential equation with an additional ${\dot{x}}^2$ term and has been used to describe self-sustained oscillations. Polyanin and Zaitsev (2002) listed some particular solutions of this equation, and Kamke (1977) gave a particular solution of this equation.

We investigated this system from mechanics and classical open systems point of view. We have considered contact Lagrangian and contact Hamiltonian formulation for one-dimensional mechanical systems described by the Levinson–Smith equation. Such a formulation enables us to obtain a geometrical description of mechanics in which thermodynamic features are naturally embedded. Hence this provides a richer and deeper understanding of the inherent features of such dynamical equations. We have illustrated the procedure with a number of examples such as the relativistic oscillator and the generalized damped oscillator. Issues related to integrability via the Chiellini condition have also been addressed and the rôle of Jacobi’s Last multiplier is illustrated. We have made use of metriplectic geometry (or GENERIC) formalism to study generalized Levinson–Smith equation which cannot be described via contact geometry.

Our next aim is to study 3D systems or cubic systems using metriplectic or GENERIC geometry in a more detailed way. Earlier we studied (Esen and Ghose 2017) this class of systems using Nambu mechanics with metriplectic structure. Many of this class of systems, for example the Rabinovich system, the Hindmarsh-Rose model etc., play an important rôle in dynamical systems. There are many open questions waiting to be answered, in fact, the GENERIC formulation of 3D dissipative systems must be studied more intensively from a more geometrical viewpoint.