SIRV fractional epidemic model of influenza with vaccine game theory and stability analysis

1.   Introduction

Most of the physical systems of interest in physics admit a description in terms of a variational principle: physical solutions are the extrema of some functional defined on the state of all possible evolutions of the system, called the action. If the action is defined in terms of the integral of a Lagrangian, then its extrema are precisely the solutions to the Euler-Lagrange equations of the Lagrangian. For mechanical systems, i.e., systems in which solutions are functions only of time, the phase space can be equipped with a symplectic structure. From this point of view, time evolution is just the flow generated by the Hamiltonian of the system, and one has at one's disposal all of the tools known from symplectic mechanics: Poisson bracket, Noether's theorem, etc. The geometric framework can be generalized to field theories by using, for instance, the multisymplectic formalism.

This description, nevertheless, excludes a large class of systems, namely, dissipative systems. In some cases, the phase space for such systems is naturally equipped with a contact structure, which can, in many ways, be seen as the odd dimensional analogue of a symplectic structure (see ^[1] for a more general discussion). What in the symplectic world were once conservation laws, now become dissipation laws ^[2]. Contact structures have made appearances in various fields in recent years, including reversible and non-reversible thermodynamics ^[3,4,5], quantum mechanics ^[6], statistical mechanics ^[7], cosmology ^[8,9] and electromagnetism ^[10]. The contact framework is well understood for mechanical systems; see ^{[2,11,12,13,14,15,16]}. The field theory analog, multicontact geometry is under current development, see ^{[17,18,19,20]} for recent efforts. One of the main challenges is a successful understanding of singular and higher-order theories.

In this work we study one such theory, namely a dissipative version of Einstein gravity. To circumvent the difficulties coming from the as of now not fully understood contact formalism, we make use of the fact that contact systems, when seen from the Lagrangian point of view, can also be formulated in terms of a variational principle, i.e., the Herglotz variational principle. The Lagrangians that fit in this framework are called action-dependent. We apply variational calculus to derive the analog of the Euler-Lagrange equations (the Herglotz equations) for this system. This is of relevance since examples of singular or higher-order dissipative field theories are scarce in the literature.

The result obtained is also relevant to the study of modifications of Einstein's theory of gravity, which would explain some observations of cosmological phenomena that do not fit within the current picture, as well as open avenues toward the successful quantization of gravity. A survey of theories of this kind is in ^[21] and ^[22]. In ^[8,23], the same Lagrangian we introduce was studied. We frame it within the broader context of the Herglotz variational principle and dissipative theories. In this sense this work is complementary to ^[8], as it clarifies the geometric nature of the objects at play and presents a set of equations that is Lorentz invariant, as opposed to the ones originally derived. This issue is also remedied in ^[23].

The work is organized as follows. In 2, we introduce the Herglotz variational principle and show how it can be equivalently formulated as a constrained optimization problem. This allows one to use the calculus of variations to derive the correct Herglotz equations of motion, which is especially relevant for field theories. In 3, we apply this language to a dissipative version of the Einstein-Hilbert Lagrangian to derive its field equations. In 4, we discuss how these equations differ from the ones originally derived in ^[8] and why they are a Lorentz-invariant generalization of them.

This article is the result of work done in ^[24].

2.   Herglotz variational problem

This chapter presents the theory of action-dependent Lagrangians. The main appeal of this formalism is that it allows for the description of non-conservative systems in terms of a variational principle, which is, in general, not possible with standard Lagrangian mechanics. The problem of finding the stationary paths of the action given by a Lagrangian of this sort is known as the Herglotz problem ^[25]. The main difficulty of this variational problem is that, as opposed to the standard variational problem of Lagrangian mechanics, it is an implicit optimization problem.

The phase space of an action-dependent Lagrangian can be equipped with a contact structure. Hence, from the Hamiltonian point of view, contact geometry is the natural framework to describe dissipative dynamics. This is well understood for mechanics, but not mature enough for field theories, and particularly for second-order theories like the Hilbert-Einstein Lagrangian. This work will focus on the Lagrangian picture and the variational methods.

There are several ways of deriving the equations of motion of the Herglotz variational principle in mechanics. The original version defines a functional on trajectories in terms of the solution to a differential equation determined by the trajectory. We refer to this as the implicit approach. Alternatively, one can implement the action dependence as a non-holonomic constraint on a standard variational problem defined on a larger configuration space and use standard variational methods. There are two distinct ways of implementing non-holonomic constraints, which are referred to as the vakonomic method and the non-holonomic method. For the Herglotz principle in mechanical systems, they are shown to be equivalent in ^[26], in the sense that they lead to the same equations.

The implicit approach to the Herglotz principle for field theories presents a difficulty because the differential equation that needs to be solved is now a partial differential equation. Nevertheless, this has been successfully done for a particular class of first-order field theories in ^[19,27]. Here, we instead follow ^[26] and use the constrained approach.

This is, to the best of the authors' knowledge, the first time that the equations of motion for a second-order action-dependent field theory have been derived. Hence, although the resulting equations are physically and geometrically sound, they cannot be compared to other results of this sort. This is relevant because, in general, the vakonomic and non-holonomic methods are not equivalent, and only one of them leads to the desired result ^[28]. This is clarified by the authors who worked in collaboration with M. Lainz and X. Rivas in ^[29]. The key result is that, when the action dependence is closed, then both methods are equivalent.

We now present the Herglotz principle for mechanical systems and first-order field theories, and show how the Herglotz equations are derived using the vakonomic method.

2.1.   Herglotz variational problem as constrained optimization

An action-dependent Lagrangian is defined on the configuration space of a non-dissipative system expanded with an extra degree of freedom. This additional degree of freedom is interpreted on-shell as the action.

In detail, consider $Q\times\mathbb{R}$, where $Q$ is the configuration space which is enlarged by an extra dimension. The Lagrangian is defined as a function $L \colon T(Q\times \mathbb{R})$ that is only zeroth-order on $z$, that is, if $(q^i, z)$ is a local chart of $Q\times \mathbb{R}$ and $(q^i, z, \dot{q}^i, \dot z)$ is the corresponding local trivialization of $T(Q\times \mathbb{R})$, then $L$ does not depend on $\dot{z}$ (or, equivalently, $dL$ annihilates the vertical vector field $\frac{\partial}{\partial \dot z}$). The constraint one imposes is

thus, for trajectories that satisfy the constraint, we have

and, indeed, $z$ tracks the action along the path, as claimed.

Let $\Omega(I, q_a, q_b, s_a)$ be the set of curves $(q, z):I = [a, b]\times\mathbb{R}$ such that $q(a) = q_a$, $q(b) = q_b$ and $z(a) = s_a$. The Herglotz problem is to determine the extrema of the functional

subject to Equation (2.1). For trajectories that satisy the constraint, we have

which resembles the classical expression of the action.

This optimization problem can be formulated equivalently by using the method of Lagrange multipliers ^[26]. Consider the Lagrangian

where $\lambda$ is the Lagrange multiplier. Then, the extrema of $S$ subject to Equation (2.1) will be the unconstrained extrema of

Because we are looking for unconstrained extrema, they will be the solutions of the Euler-Lagrange equations for $\tilde{L}$. The equation for $z$ is

or, equivalently,

The Euler-Lagrange equations for the other coordinates are

and, after incorporating Equation (2.5) and dividing through by $\lambda$, one obtains

These are the Herglotz equations.

2.2.   Action-dependent field theory

We now introduce the Herglotz problem for field theories and derive the corresponding Herglotz equations.

2.2.1.   Classical field theory and Lagrangian densities

The passage from mechanics to field theory requires some care. The given data are usually some smooth fiber bundle $E$ over a base $M$ of dimension $n$, which we assume to be orientable and hence endowed with at least a volume form. The base $M$ is usually, but not always, taken to represent spacetime. Field configurations are sections of this bundle, and the values that the field takes are modeled by the fiber of $E$. The basic problem is to identify the configurations that are extrema of a given action functional $S \colon \Gamma(E) \to \mathbb{R}$. The action is usually written as the integral of a Lagrangian over a region $D \subseteq M$ of the base, i.e., a Lagrangian is some sort of map from field configurations to top forms of $M$, i.e., $\mathcal{L} \colon \Gamma(E) \to \Omega^{n}(M)$, such that

One of the fundamental constraints on $\mathcal{L}$ is that it must be local, i.e., $\mathcal{L}(\phi)_p$ should only depend on the value of the field $\phi$ at $p$ and a finite number of its derivatives at $p$. In other words, $\mathcal{L}$ is to be a bundle map from the $k$-th jet bundle of $E$, $J^kE$, to the bundle of top forms $\bigwedge^{n}T^\ast M$ such that if $j^k \phi \in \Gamma(J^kE)$ is the prolongation of some field configuration $\phi \in \Gamma(E)$, then

The integer $k$ is called the order of the Lagrangian. Volume forms are one-dimensional, which means that, for a given choice of coordinates of the base, $x^\mu$, then there exists a unique $L\colon J^1E \to \mathbb{R}$ such that, on the coordinate domain,

where $\, {\mathrm{d}}^nx$ is the local volume form of $M$ induced by the coordinates.

Using the calculus of variations, one can show that the stationary configurations of an action functional defined by a first-order Lagrangian satisfy the Euler-Lagrange equations of field theory:

This expression makes sense given the choice of a local trivialization of $E$, i.e., $(x^\mu, \phi^a)$, which gives rise to a local trivialization of $J^1E$, i.e., $(x^\mu, \phi^a, \phi^a_\mu)$. Note that the Einstein summation convention is assumed from this point on, unless otherwise stated.

2.2.2.   Action flux

We now wish to generalize this description to account for action-dependent Lagrangians. A cursory look at the Herglotz equations would suggest a field theory analog of the form

These equations have also been proposed in the literature ^[18,20,27]. The question then becomes what should be the geometric nature of $z$. In the language of bundles that we have just introduced, the constraint in Equation (2.1) becomes

The field theory version should then be

where $D\subset M$ is an $n$-dimensional submanifold of $M$ over which we wish to extremize the action. It is now clear that $z$ must be a form of degree $n-1$ so that it can be integrated over submanifolds of the base $M$ of codimension 1. In other words, $z$ is the action flux. The differential version of Equation (2.11) is then

This is analogous to Equation (2.1).

There is, by way of contraction with a volume form, an isomorphism between $(n-1)$-forms and vector fields such that the exterior derivative becomes the divergence. In particular, a choice of coordinates $x^\mu$ on the base induces a trivialization $(x^\mu, z^\nu)$ on $\bigwedge^{n-1} T^\ast M$, such that, for $\alpha \in \bigwedge T^\ast _p M$,

where the symbol $\lrcorner$ denotes the contraction of a tangent vector with a form. Then, for a differential form $z$ of degree $n-1$ whose components in these coordinates are $z^\nu$, it holds that

This means that the coordinate expression of Equation (2.13) is

For some bundle $E \to M$, the corresponding Herglotz problem is formulated in the enlarged bundle $E \oplus \bigwedge^{n-1}T^\ast M \to M$, where $\oplus$ is the Whitney sum. Consider a Lagrangian of the form $\mathcal{L} \colon J^kE \oplus \bigwedge^{n-1}T^\ast M \to \bigwedge^n T^\ast M$ (so that, crucially, $\mathcal{L}$ does not depend on any of the derivatives of the action flux). The Herglotz problem for field theory is then to find the sections $(\phi, z)$ that extremize the functional $S(\phi, z) = \int_{\partial D} z$ which is subject to the constraint $\, {\mathrm{d}} z = \mathcal{L}$. If $(\phi, z)$ is one such section, then

and we can interpret $S$ as the action.

2.2.3.   Constrained optimization in field theory

Just like before, we turn this constrained optimization problem into an unconstrained one by using Lagrange multipliers. The expanded action for a first-order Lagrangian, similar to Equation (2.6), is

Let us write down the integrand of Equation (2.17) as an expanded Lagrangian:

Note that $\tilde{ \mathcal{L}}$ is now the Lagrangian for a theory defined on the expanded bundle $J^1E \oplus J^1\bigwedge^{n-1}T^\ast M \to M$, so $z$ is a dynamical degree of freedom.

Given that the Lagrangian is of first order, extrema of this action functional will be solutions to the Euler-Lagrange equations for this Lagrangian, which become the Herglotz equations upon imposing the constraint. We present them in the next section. Nevertheless, there is nothing preventing one from calculating the explicit variation of the action, which leads to the equations of motion for a theory of any order. This is the approach we follow in the next chapter.

2.2.4.   Herglotz equations for field theory

Finally, we derive the Herglotz equations for field theory from the expanded Lagrangian in Equation (2.18). The equations for the action flux are

where $(x_\mu, z^\nu, z^\nu_\mu)$ is the trivialization of $J^1(\bigwedge^{n-1}T^\ast M)$ induced by the choice of coordinates on the base, as defined by Equation (2.14). Rearranging, one obtains

This equation actually constrains the type of action dependence that is allowed in $L$. We will see later that, in the context of relativity, it forces the dissipation form to be closed.

The equations for the field are

and, after substituting in Equation (2.19) and dividing through by $\lambda$, we arrive at the field theory Herglotz equations

3.   Action-dependent Einstein gravity

In this chapter, we apply the language and tools developed in the previous chapter to the specific case of Einstein gravity. We first describe the Lagrangian from which the Einstein field equations arise. Then, introduce an action-dependent version of it and we derive its field equations.

3.1.   Einstein-Hilbert Lagrangian

As is well known, the Einstein field equations can be obtained from a variational principle. The classical Lagrangian that gives rise to these equations is the Einstein-Hilbert Lagrangian. We formulate it in the language of Section 2.2.1. The field of interest in relativity is the metric on the given spacetime $M$, which we take to be of dimension 4 and orientable. Hence, the bundle of interest to us is a subbundle of the second symmetric power of the cotangent bundle of $M$,

Specifically, it is the subbundle determined by the condition of non-degeneracy. We denote it by $G(M) \to M$.

As advertised, the theory of general relativity is a second-order theory, which means that the Einstein-Hilbert Lagrangian must be a bundle map from $J^2 G(M)$ to $\bigwedge^n T^\ast M$. Specifically, given a metric $g \in \Gamma(G(M))$,

Here, we use $\omega_ g$ to denote the volume form determined by $g$, which, in a choice of coordinates $x^\mu$, becomes

where $\sqrt{g}$ is the square root of the absolute value of the determinant of the expression of the metric in the coordinates $x^\mu$. The other factor, $R(g)$, is the scalar curvature of $g$, which is defined as the trace of the Ricci tensor:

The Ricci tensor is of type $(0, 2)$; so by contracting with $g^{-1}$, i.e., the metric induced on $T^\ast M$, we obtain a $(1, 1)$ tensor, whose trace is well defined. The coordinate expression of the components of the Ricci tensor, $R_{ab}$, is

where $\Gamma_{ab}^c$ are the Christoffel symbols of the Levi-Civita connection determined by $g$. These contain the first derivatives of the metric, so $R_{ab}$ and, hence, $R$ contain the second derivatives of the metric, and the Einstein-Hilbert Lagrangian is indeed of second order.

The Einstein-Hilbert action is therefore

for some domain $D$ on which the integral is finite. A variation of this action leads one to the Einstein field equations

More precisely, these are the Einstein field equations in a vacuum, since one can add various matter terms to the Einstein-Hilbert Lagrangian which lead to the Einstein equations in the presence of matter;

The object $T_{ab}$ is the energy-momentum tensor, and it collects all of the terms coming from the presence of matter. See §4 of ^[30] for a detailed derivation.

3.2.   An action dependent Einstein-Hilbert Lagrangian

What kind of action dependence can we incorporate into the Einstein-Hilbert Lagrangian? The simplest one is a linear dissipation term:

Now, $\mathcal{L}_{\text{E-H}}$ is defined on the expanded bundle, $G(M) \oplus \bigwedge^3T^\ast M \to M$. Hence, $\theta$ must be a 1-form on $M$, which we will refer to as the dissipation form.

The coordinate expression of this dissipation term is

where, once again, we use the coordinates defined in Equation (2.14). Then, Equation (3.9) becomes

This Lagrangian does not exactly match the one proposed in Equation (9) of ^[8]. The discrepancy is down to a different choice of coordinates. Indeed, in the previous computation, we used the isomorphism between $\Omega^3(M)$ and $\Gamma(TM)$ induced by contracting with ${\mathrm{d}}^4x$. Instead, we may contract with the volume form induced by the metric, $\omega_g$. Let $\zeta^\mu$ be the components of $z$ in this new choice of coordinates, i.e.,

which implies that $z^\mu = \sqrt{g} \zeta^\mu$. Given these new coordinates, Equation (3.10) looks like

This is the Lagrangian proposed in Equation (9) of ^[8].

We now write down the constraint in Equation (2.13) for this Lagrangian. In the original coordinates for the action flux, we have

so

In the other set of coordinates, induced by contracting with $\omega_g$, one sees

where $\nabla$ is the covariant derivative induced by $g$. We have made use of a useful identity about the divergence:

This is the statement that the divergence induced by the volume form of a metric coincides with the trace of the covariant derivative.

Given the new coordinates, the constraint becomes

which is the same form that appears in Equation (8) of ^[8].

3.3.   Derivation of the field equations

We now apply the method of Lagrange multipliers, as described in the previous chapter, to derive a modified version of Einstein's equations. The expanded Lagrangian is

We will compute the variation of the corresponding expanded action, $\tilde{S}(g, z) = \int_D \tilde{L}_\text{E-H} \circ (j^2g, j^1z)$, with respect to the two dynamical degrees of freedom, $z$ and $g$.

3.3.1.   Variation of the action flux

The variation with respect to the action flux is

The first integral is a boundary term coming from an integration by parts. It vanishes if we assume that the variations vanish at the boundary of $D$. If the action is stationary, then its variation must vanish for any variation of the fields. This means that the second term of Equation (3.15) must vanish, since, in particular, we may choose not to vary the metric. Hence, the quantity inside the brackets must vanish since it vanishes when integrated against any variation. Therefore,

In other words, $\, {\mathrm{d}} \lambda = \lambda \theta$. As we had advertised before, this forces the dissipation form $\theta$ to be closed, because

hence,

and we conclude that $\, {\mathrm{d}} \theta = 0$ provided that $\lambda$ does not vanish.

3.3.2.   Variation of the metric

We continue the calculation from Equation (3.15). We may now only consider the last integral, as we can vary $g$ and $z$ independently. We will follow the derivation in ^[30] for as long as we can. In particular, we take the spacetime $M$ to be closed, and hence avoid consideration of Gibbons-Hawking-York-type boundary terms. Since $R\sqrt{g} = g^{ab}R_{ab}\sqrt{g}$, from the product rule, its variation results in three terms as follows:

The first term is already in the form required to apply the fundamental theorem of the calculus of variations. The third one uses the standard result:

The first and third terms of Equation (3.17) can be combined into

In the standard derivation of Einstein's equations, one shows that the middle integral of Equation (3.17) actually vanishes; thus, if Equation (3.18) is to vanish for any variation $\delta g_{ab}$, or, equivalently, for any variation of the inverse metric $\delta g^{ab}$, the integrand of Equation (3.18) itself must vanish. This gives Einstein's equations. In the presence of $\lambda$, however, the middle integral does not vanish and it contributes additional terms to the equations.

We compute the variation of the middle integral in Equation (3.17). The variation of the Ricci curvature can be shown to be

so

and, if $\lambda$ were not there, this integral would vanish because of the divergence theorem and the fact that the variations vanish on the boundary of $D$. In the presence of $\lambda$, we perform an integration by parts:

The first integral vanishes because it is the integral of a divergence and the variations vanish on the boundary of $D$. The second integral is the source of the additional terms. We split it into two terms.

The variation of the Christoffel symbols can be shown to be

Using this and Equation (3.16) (since $\nabla_n \lambda = \partial_n \lambda$), we compute the following for the first integral:

The presence of $g^{ab}$ means that the indices $a$ and $b$ are symmetrized, so

This means that Equation (3.21) simplifies to

Let us perform an integration by parts for the first integral. Introducing the shorthand $X^{bn} = g^{ab}g^{nk}\delta g_{ak}$, we compute

so

The first integral is the integral of a divergence, so it vanishes. We are left with the second which we can expand into

As a last step, we use the identity

to write our integral as a variation with respect to the inverse metric.

Without going through the details again, the other integral in Equation (3.22) can be brought to the form

There is still another integral we need to evaluate, and it is the second term in the variation of $R_{ab}$, namely

Because $m$ and $k$ are symmetrised, the second and third terms cancel, leaving us with

We have calculated all of the integrals that we need. Before we put them all together, let us make the following observation:

which uses the fact that $\theta$ must be closed. We may therefore define the following (0, 2) symmetric tensor:

whose components are

where parentheses surrounding indices indicate symmetrization. All three expressions are equal because $\nabla_a \theta_b = \nabla_b \theta_a$. Nevertheless, we will use the second one to make the symmetry of the indices explicit. So, after liberal relabeling of the indices, we find that Equation (3.15) becomes

with $K_{ab}$ defined as in Equation (3.27) and $K = g^{mn}K_{mn}$ as its trace.

Applying the fundamental theorem of the calculus of variations, the action will be stationary if and only if the integrands of both terms vanish. From the first integral, we get Equation (3.16), which we have already used. And, from the second one, we get the modified Einstein field equations

These equations coincide with the ones derived in ^[23].

4.   Significance of the equations

In this chapter, we discuss the equations that we have obtained and how they compare to those appearing in existing publications. We also make the case that the version we have derived is a more adequate version.

4.1.   Dissipation tensor

Let us recap what we did in the previous chapter. We have shown, by computing the variation of the corresponding action, that the field equations of an Einstein-Hilbert Lagrangian with linear dissipation, namely,

are

where $K_{ab}$ represents the components of the $(0, 2)$ symmetric tensor

We will call $K$ the dissipation tensor. Since the first two terms of Equation (4.2) have zero divergence (i.e., they are the components of the Einstein tensor), it must be the case that, on-shell, the divergence of the second two terms also vanishes. This imposes a constraint on the space of solutions to Equation (3.30), which depends on the dissipation form $\theta$. Namely,

This means that, if we were to couple a matter term to Equation (3.9), its energy-momentum tensor need not, in general, have zero divergence. Specifically, what must have zero divergence will be a combination of the energy-momentum tensors of the matter fields and terms containing the dissipation 1-form. Nevertheless, further investigation is required to determine the precise way in which the dissipation tensor governs the non-conservation of other quantities.

4.2.   Non-covariance of existing equations

These equations are not the ones obtained in ^[8]. For the same Lagrangian, the equations derived are

where $\tilde{K} = g^{ab}\tilde{K}_{ab}$ and $\tilde{K}_{ab}$ is

These cannot possibly represent the components of a tensor. Very explicitly, for the flat Minkowski metric, their expression in Cartesian coordinates is 0. If they represented the components of a tensor, then they would also vanish for any other choice of coordinates for the flat metric. Nevertheless, in spherical coordinates, one computes

This means that, for example,

which is certainly non-zero if $\theta_t$ does not vanish. Hence, the object derived in ^[8] is not coordinate-independent, so it cannot possibly represent meaningful physics.

There is another fact that points to the equations in ^[8] not being what one would expect as the Herglotz equations coming from a second-order action-dependent Lagrangian. The Herglotz equations for the harmonic oscillator with linear dissipation lead to equations that are linear in the dissipation coefficient (^[2]). However, the Lagrangian for this system is first-order, whereas, as we had already discussed, the Einstein-Hilbert Lagrangian is actually second-order. There is a second-order Lagrangian with linear dissipation, called the damped Pais-Uhlenbeck oscillator, whose equations of motion are derived in ^[15]. These are, in fact, not linear in the dissipation coefficient, but, rather, quadratic. In our case, the dissipation form plays the role of the dissipation coefficient and, indeed, ${\bf{K}}$ is quadratic in it. The equations in ^[8] instead lack a quadratic term.

One can pinpoint the exact reason for the problems with Equation (4.5). One of the simplifying assumptions made in their derivation is to only consider certain terms of the Ricci curvature. Specifically, the Ricci curvature consists of four terms. Two of them are contractions of the Christoffel symbols with themselves, and the other two are derivatives of the Christoffel symbols. In the classical case, without dissipation, one can show that the second two terms are actually a divergence, so they do not contribute to the variation of the Einstein-Hilbert action; and the resulting equations remain unchanged (see ^[31,32]). For an action-dependent theory, however, adding a divergence to the unexpanded Lagrangian does not lead, in general, to the same equations ^[23,33].

5.   Conclusions

The are three main ideas presented in this article.

First, we showed how the Herglotz problem can be turned from a constrained optimization problem to an unconstrained one by promoting the action dependence to a dynamic degree of freedom and using Lagrange multipliers to implement the non-holonomic constraint.

Second, we described how the Einstein-Hilbert Lagrangian can be modified with an action dependence in a coordinate-independent manner. This allows one to derive a correct, Lorentz-invariant set of field equations that remedy the issues present in previous derivations.

Finally, the computation performed constitutes an important example for the ongoing development of contact geometry and its applications, since it is a singular second-order field theory. Having a concrete example at hand will aid in understanding these systems.

There are various avenues for future follow-up work. One can consider more general dissipation terms to add to the Einstein-Hilbert Lagrangian, as well as study their phenomenology. It will also be interesting to consider the boundary effects in the case of manifolds with a boundary (the appropriate Gibbons-Hawking-York term).

General relativity has several equivalent formulations ^[34,35]. It would be interesting to add dissipation to these formalisms and study their properties and relations.

Finally, the current tools in contact geometry fall short of completely describing this kind of Lagrangians. A more general geometric structure, akin to multisymplectic geometry, needs to be developed for more general action-dependent Lagrangians in order to describe relevant theories. In this line, the multicontact structure recently presented in ^[20] could be the adequate geometric framework for action-dependent gravity.

Acknowledgments

The authors wish to thank the reviewers for their constructive suggestions which improved the final version of the paper.

Jordi Gaset also acknowledges financial support from the Ministerio de Economía y Competitividad. (Spain), project MTM2014-54855-P, and from the Ministerio de Ciencia, Innovación y Universidades (Spain), projects PGC2018-098265-B-C33 and D2021-125515NB-21.