Research article

Simulated biomechanical effect of aspheric transition zone ablation profiles after conventional hyperopia refractive surgery

  • Received: 06 January 2021 Accepted: 08 March 2021 Published: 12 March 2021
  • We studied the effects of the aspheric transition zone on the optical wavefront aberrations, corneal surface displacement, and stress induced by the biomechanical properties of the cornea after conventional laser in situ keratomileusis (LASIK) refractive surgery. The findings in this study can help improve visual quality after refractive surgery. Hyperopia correction in 1-5D was simulated using five types of aspheric transition zones with finite element modeling. The algorithm for the simulations was designed according to the optical path difference. Wavefront aberrations were calculated from the displacements on the anterior and posterior corneal surfaces. The vertex displacements and stress on the corneal surface were also evaluated. The results showed that the aspheric transition zone has an effect on the postoperative visual quality. The main wavefront aberrations on the anterior corneal surface are defocus, y-primary astigmatism, x-coma, and spherical aberrations. The wavefront aberrations on the corneal posterior surface were relatively small and vertex displacements on the posterior corneal surface were not significantly affected by the aspheric transition zone. Stress analysis revealed that the stress on the cutting edge of the anterior corneal surface decreased with the number of aspheric transition zone increased, and profile #1 resulted in the maximum stress. The stress on the posterior surface of the cornea was more concentrated in the central region and was less than that on the anterior corneal surface overall. The results showed that the aspheric transition zone has an effect on postoperative aberrations, but wavefront aberrations cannot be eliminated. In addition, the aspheric transition zone influences the postoperative biomechanical properties of the cornea, which significantly affect the postoperative visual quality.

    Citation: Ruirui Du, Lihua Fang, Binhui Guo, Yinyu Song, Huirong Xiao, Xinliang Xu, Xingdao He. Simulated biomechanical effect of aspheric transition zone ablation profiles after conventional hyperopia refractive surgery[J]. Mathematical Biosciences and Engineering, 2021, 18(3): 2442-2454. doi: 10.3934/mbe.2021124

    Related Papers:

    [1] Elias Maciel, Inocencio Ortiz, Christian E. Schaerer . A Herglotz-based integrator for nonholonomic mechanical systems. Journal of Geometric Mechanics, 2023, 15(1): 287-318. doi: 10.3934/jgm.2023012
    [2] Xavier Rivas, Daniel Torres . Lagrangian–Hamiltonian formalism for cocontact systems. Journal of Geometric Mechanics, 2023, 15(1): 1-26. doi: 10.3934/jgm.2023001
    [3] Jordi Gaset, Arnau Mas . A variational derivation of the field equations of an action-dependent Einstein-Hilbert Lagrangian. Journal of Geometric Mechanics, 2023, 15(1): 357-374. doi: 10.3934/jgm.2023014
    [4] Alexander Pigazzini, Cenap Özel, Saeid Jafari, Richard Pincak, Andrew DeBenedictis . A family of special case sequential warped-product manifolds. Journal of Geometric Mechanics, 2023, 15(1): 116-127. doi: 10.3934/jgm.2023006
    [5] Marcin Zając . The dressing field method in gauge theories - geometric approach. Journal of Geometric Mechanics, 2023, 15(1): 128-146. doi: 10.3934/jgm.2023007
    [6] Robert I McLachlan, Christian Offen . Backward error analysis for conjugate symplectic methods. Journal of Geometric Mechanics, 2023, 15(1): 98-115. doi: 10.3934/jgm.2023005
    [7] Álvaro Rodríguez Abella, Melvin Leok . Discrete Dirac reduction of implicit Lagrangian systems with abelian symmetry groups. Journal of Geometric Mechanics, 2023, 15(1): 319-356. doi: 10.3934/jgm.2023013
    [8] Markus Schlarb . A multi-parameter family of metrics on stiefel manifolds and applications. Journal of Geometric Mechanics, 2023, 15(1): 147-187. doi: 10.3934/jgm.2023008
    [9] Valentin Duruisseaux, Melvin Leok . Time-adaptive Lagrangian variational integrators for accelerated optimization. Journal of Geometric Mechanics, 2023, 15(1): 224-255. doi: 10.3934/jgm.2023010
    [10] Francesco Bonechi, Jian Qiu, Marco Tarlini . Generalised Kähler structure on CP2 and elliptic functions. Journal of Geometric Mechanics, 2023, 15(1): 188-223. doi: 10.3934/jgm.2023009
  • We studied the effects of the aspheric transition zone on the optical wavefront aberrations, corneal surface displacement, and stress induced by the biomechanical properties of the cornea after conventional laser in situ keratomileusis (LASIK) refractive surgery. The findings in this study can help improve visual quality after refractive surgery. Hyperopia correction in 1-5D was simulated using five types of aspheric transition zones with finite element modeling. The algorithm for the simulations was designed according to the optical path difference. Wavefront aberrations were calculated from the displacements on the anterior and posterior corneal surfaces. The vertex displacements and stress on the corneal surface were also evaluated. The results showed that the aspheric transition zone has an effect on the postoperative visual quality. The main wavefront aberrations on the anterior corneal surface are defocus, y-primary astigmatism, x-coma, and spherical aberrations. The wavefront aberrations on the corneal posterior surface were relatively small and vertex displacements on the posterior corneal surface were not significantly affected by the aspheric transition zone. Stress analysis revealed that the stress on the cutting edge of the anterior corneal surface decreased with the number of aspheric transition zone increased, and profile #1 resulted in the maximum stress. The stress on the posterior surface of the cornea was more concentrated in the central region and was less than that on the anterior corneal surface overall. The results showed that the aspheric transition zone has an effect on postoperative aberrations, but wavefront aberrations cannot be eliminated. In addition, the aspheric transition zone influences the postoperative biomechanical properties of the cornea, which significantly affect the postoperative visual quality.



    An invariant volume is a powerful tool for understanding the asymptotic nature of a dynamical system. In particular, it is a well-known fact that unconstrained mechanical systems are volume-preserving. The case of nonholonomic systems is more nuanced as these systems generally fail to preserve the symplectic form (which follows from the fact that nonholonomic systems are not variational) and hence, the induced volume form. This makes the study of invariant volumes in nonholonomic systems nontrivial and the existence of an invariant volume is an exceptional scenario which should not be expected in general [1,2]. A famous example of this is the Chaplygin sleigh; this system, although energy-preserving, experiences volume "dissipation" which results in asymptotic stability (cf. [3] for a general discussion on stability of nonholonomic systems or [4] for an interpretation via impact systems).

    The existence of an invariant volume for a nonholonomic system offers two key insights. The first is the usual case in dynamical systems where an invariant volume allows for the use of the Birkhoff Ergodic Theorem (cf. e.g. 4.1.2 in [5]), for recurrence (with the caveat of the volume being finite although some results still hold for infinite measures [6]), and as an obstruction to asymptotic stability. The other is unique to nonholonomic systems; even though nonholonomic systems are not Hamiltonian, "nonholonomic systems which do preserve volume are in a quantifiable sense closer to Hamiltonian systems than their volume changing counterparts, " [7] (see also [8,9,10,11,12]). Therefore, being able to find an invariant volume for a nonholonomic system allows for understanding its asymptotic behavior and can provide a way to "Hamiltonize" a nonholonomic system (although most systems with an invariant volume still cannot be Hamiltonized, cf. [13]).

    There exists an abundance of research into finding invariant volumes for nonholonomic systems where the constraints are linear and symmetries are present: Chaplygin systems are studied in, e.g. [11,14,15,16,17,18,19], Euler-Poincaré-Suslov systems are studied in, e.g. [20,21], systems with internal degrees of freedom are studied in, e.g. [20,22,23], and [24] studies the case of symmetric kinetic systems where the dimension assumption does not hold. Related work on asympotic dynamics may be found in [25]. There are also results demonstrating that no invariant volumes exist, e.g. [26,27].

    The contribution of this work is to consider conditions for the existence of invariant volume in nonholonomic systems such that

    1. the constraints need not be linear/affine in the velocities,

    2. the analysis can be carried out on the ambient manifold rather than resorting to local coordinates, and

    3. absolutely no symmetry assumptions are used.

    Such an approach seems to be new.

    Our main result is an existence condition for an invariant measure for a hyper-regular Lagrangian system with an admissible set of nonlinear nonholonomic constraints obeying Chetaev's rule. In the linear case the existence of a basic invariant density reduces to checking whether or not a certain 1-form can be made to be exact, [28], and this special case is similar to results in [15,24].

    Preliminaries on nonholonomic systems are presented in Section 2. Section 3 presents the construction of the "extended nonholonomic vector fields" on the entire cotangent space TQ which restricts to the nonholonomic vector field on M. Although there does not exist a canonical choice for the extended vector field, they will agree once restricted. An advantage of this approach is that one is able to design extended vector fields which ensure that the constraint manifold is an exponentially stable invariant manifold. The divergence calculation for a nonholonomic system is performed in Section 4 via extended nonholonomic vector fields. The main results, Theorem 5.2 and Theorem 5.6, are proved in Section 5 (cf. Theorems 5.2 and 5.6). Additionally, it is shown that the invariant volume does not depend on the choice of extended nonholonomic vector field chosen and is unique up to a constant of motion. Section 6 presents an interpretation of the linear constraint case to the torsion of the nonholonomic connection, which seems to be a novel observation. This paper concludes with examples in Section 7.

    This paper is a continuation of the work done in [29] and, as such, many of the results below can be found there.

    We will first briefly cover the case of unconstrained mechanical systems before discussing nonholonomic systems. A smooth (finite-dimensional, dim(Q)=n) manifold Q is called the configuration space, the tangent bundle TQ is called the state space, and the cotangent bundle TQ is called the phase space. The bundle projection maps will be denoted by τQ:TQQ and πQ:TQQ.

    Lagrangian systems are described by a smooth real-valued function on the state space called the Lagrangian, L:TQQ. The dynamics generated from the Lagrangian function arise from Hamilton's principle and are given by the Euler-Lagrange equations,

    ddtL˙qLq=0.

    In constrast to Lagrangian systems, Hamiltonian systems evolve on the phase space. For a given Lagrangian, denote FL:TQTQ as its fiber derivative.

    Definition 2.1. A Lagrangian, L:TQR, is hyperregular if the fiber derivative is a global diffeomorphism. Furthermore, a Lagrangian is natural if it has the form

    L(v)=12g(v,v)τQV(v),

    where g is a Riemannian metric on Q and VC(Q).

    Throughout this work, all Lagrangians will be assumed to be hyperregular. Moreover, when affine constraints are considered, the Lagrangian will be assumed to be natural.

    For a given Lagrangian, we define the Hamiltonian H:TQR via the Legendre transform:

    H(p)=p,vL(v),p=FL(v).

    Let ω=dqidpi be the canonical symplectic form on TQ. The resulting Hamiltonian vector field arising from the Hamiltonian H and the symplectic form ω is denoted by XH and is given by

    iXHω=dH,

    where iXω=ω(X,) is the contraction. An important feature of Hamiltonian vector fields is that they preserve the symplectic form, and thus are volume-preserving.

    Theorem 2.2 (Liouville's Theorem). Hamiltonian dynamics preserve the symplectic form and, additionally, preserve the volume form ωn.

    Liouville's theorem does not apply to general nonholonomic systems as they are generally not symplectic. Determining when an analogous version of this theorem applies to constrained systems is the main goal of this work.

    For unconstrained Lagrangian systems, motion is allowed to occur in the entire state space, TQ. For constrained systems, motion is restricted to lie on a constraint manifold NTQ and we tacitly assume that τQ(N)=Q. For the purposes of this work, we will not differentiate between nonholonomic and holonomic constraints as the machinery for the former will suffice for the latter.

    Let NTQ be the constraint submanifold which is, locally, described via the zero level-set of a regular function G:TQRk, k<n=dim(Q). Let us denote each component of this function by Ψα, α=1,,k, i.e.

    N=kα=1(Ψα)1(0).

    To transfer to the Hamiltonian point of view, we assume that the Lagrangian, L:TQR, is hyper-regular. The constraints on the cotangent bundle become

    Φα=(FL1)Ψα,M=FL(N)=kα=1(Φα)1(0).

    These constraints are admissible if the following matrix is non-singular [18].

    mαβ=CdΦβ(XΦα),(mαβ)=(mαβ)1,

    where XΦα is the Hamiltonian vector field generated by Φα and C:T(TQ)T(TQ) is the FL-related almost-tangent structure, explained below.

    The constrained Euler-Lagrange equations become modified via Chetaev's rule (which is not necessarily the physically correct procedure, cf. [30] for a discussion), which will provide equivalent results to those in [31] where the "almost-tangent" structure of the tangent bundle is utilized. For Lagrangian systems, Chetaev's rule states that if we have the nonlinear constraints Ψα, then the constraint forces have the following form

    ddt(L˙q)Lq=λαSdΨα, (2.1)

    where S:T(TQ)T(TQ) is a (1, 1)-tensor called the almost-tangent structure and, in local coordinates, is given by

    S=˙qidqi.

    If, rather than being general nonlinear, the constraints are affine in the velocities,

    Ψα(v)=ηα(v)+τQξα(v),

    where ηαΩ1(Q) are 1-forms and ξαC(Q) are functions, then (2.1) becomes

    ddt(L˙q)Lq=λατQηα.

    However, this work will focus on the Hamiltonian formalism. The constraint manifold, MTQ, is locally described by the joint level-set of the functions Φα:TQR where Φ(p)=ΨFL1(p). Moving (2.1) to the cotangent side, the constrained Hamiltonian equations of motion become

    iXMHω|M=dH|M+λαCdΦα|M, (2.2)

    where C:T(TQ)T(TQ) is the FL-related almost-tangent structure, i.e. if L is natural, then

    C(αidqi+βjdpj)=gijβjdqi.

    Throughout this work, L:TQR will be a hyperregular Lagrangian with corresponding Hamiltonian H:TQR. The constraint manifold will be called NTQ or M=FL(N)TQ. These submanifolds will be locally described by the joint zero level-set of a collection of smooth functions,

    N=kα=1(Ψα)1(0),M=kα=1(Φα)1(0).

    If the constraints are affine in the velocities/momentum, the Lagrangian will be assumed to be natural with Riemannian metric g. In this case, the constraining functions will have the form

    Ψα(v)=ηα(v)+τQξα(v),Φα(p)=P(Wα)(p)+πQξα(p),

    where ηα are 1-forms, ξα are functions, Wα=FL1(ηα) are vector fields, and P(Wα) is the vector field's momentum

    P(Wα)(p)=p,Wα(πQ(p)).

    In this case of affine constraints, DTQ will be the distribution

    Dq=kα=1kerηαq,

    and D=FL(D)TQ.

    Finally, the nonholonomic vector fields with Hamiltonian H and constraint manifold M will be denoted by XMH.

    Given a constraint manifold, MTQ, we can determine the nonholonomic vector field, XMHX(M) via (2.2). Commonly local, noncanonical, coordinates are chosen for M (cf. §5.8 in [20] and [32]). However, we will instead work with a tubular neighborhood (for ease, this will be written as the entire manifold TQ) and define an extended vector field XextHX(TQ) such that XextH|M=XMH. This section outlines an intrinsic (albeit non-unique) way to determine such a vector field.

    Definition 3.1. For a given constraint submanifold MTQ, a realization of M is an ordered collection of functions C:={Φα:TQR} such that zero is a regular value of G=Φ1××Φk and

    M=α(Φα)1(0).

    If the functions Φα are affine in momenta, i.e. Φα=P(Xα)+πQfα, then the realization is called affine.

    Remark 3.2. In the case where the Lagrangian is natural (which provides a Riemannian metric on Q) and the constraint submanifold is affine, we can choose the realization to be affine:

    C={P(W1)+πQξ1,,P(Wk)+πQξk},

    where Wi=FL1ηi=(ηi).

    By replacing M with a realization C, we can extend the nonholonomic vector field to a vector field on TQ such that Φα are first integrals. Recall that the form of the nonholonomic vector field is iXMHω=dH+λαCdΦα. We construct the extended nonholonomic vector field, ΞCH, by requiring that:

    (NH.1) iΞCHω=dH+λαCdΦα for smooth functions λα:TQR, and

    (NH.2) LΞCHΦα=0 for all ΦαC.

    Under reasonable compatibility assumptions on C (cf. §3.4.1 in [18], see Definition 3.4 below), such a vector field exists and is unique. However, given two different realizations, C and C, of the same constraint manifold M, it is not generally true that ΞCH=ΞCH, however ΞCH|M=ΞCH|M.

    Remark 3.3. The constraint manifold is given by the joint zero level-sets of the Φα while the realization provides additional irrelevant information off of the constraint manifold. This is why ΞCHΞCH but they agree once restricted as will be proved in Proposition 3.6 below.

    Definition 3.4. For a realization C={Φ1,,Φk}, the constraint mass matrix, (mαβ), is given by "orthogonally" pairing the constraints,

    mαβ=CdΦα(XΦβ).

    The realization is admissible if this matrix is non-singular on a tubular neighborhood of M. The inverse matrix will be denoted by (mαβ)=(mαβ)1.

    When the constraints are affine (along with a natural Lagrangian), the constraint mass matrix becomes

    mαβ=g(Wα,Wβ)=ηα(Wβ),

    and admissibility of the constraints is equivalent to the constraints being linearly independent.

    We can now write down a formula for ΞCH. Using (NH.1) and (NH.2), we get that (where {,} is the standard Poisson bracket)

    LΞCHΦβ=iXΦβω(ΞCH)=iΞCHω(XΦβ)=dH(XΦβ)λαCdΦα(XΦβ)={Φβ,H}λαCdΦα(XΦβ)=0.

    This implies that {Φβ,H}=mαβλα. Due to the constraint mass matrix being nondegenerate, the multipliers are uniquely determined and the extended nonholonomic vector field is determined by

    iΞCHω=dHmαβ{H,Φα}CdΦβ (3.1)

    Definition 3.5. The 1-form on TQ given by

    νCH:=dHmαβ{H,Φα}CdΦβ,

    is called the nonholonomic 1-form with respect to the realization C.

    Proposition 3.6. Given two different realization, C and ˜C, the extended nonholonomic vector fields given by (3.1) agree on M.

    Proof. Let C={Φα} and ˜C={˜Φα} be two different realizations of the same manifold M and let xM. Their differentials span the annihilator,

    span{dΦαx}=span{d˜Φαx}=Ann(TxM).

    As the map C:T(TQ)T(TQ) is linear, we have

    span{CdΦαx}=span{Cd˜Φαx}.

    As the multipliers are unique, we must have

    mαβ{H,Φα}CdΦβ|M=˜mαβ{H,˜Φα}Cd˜Φβ|M,

    and therefore the resulting extended nonholonomic vector fields must agree on M.

    Remark 3.7. A consequence of Proposition 3.6 is that this procedure of extending the nonholonomic vector field is still valid when M is not globally defined as the level-set of functions. Suppose that U,VM are two neighborhoods characterized by the zero level-sets of functions Φ,ϕ:TQR. Then Φ and ϕ both define different realizations on UV. However, their corresponding vector fields agree on the intersection. This observation will make the constructions required for Theorem 5.2 well-defined even if M cannot be globally defined as a level set.

    The extended nonholonomic vector field, XMH, is a vector field on the ambient space TQ and has M as an invariant manifold. As such, an integral curve of XMH is only a valid trajectory if its initial condition lies in M. In practice using numerical techniques, round off errors can result in trajectories drifting away from M. To counter this effect, the condition (NH.2) can be replaced to include something in the spirit of a "control Lyapunov function" [33]. This modification will only be explored in this section while the rest of this paper will be concerned with the unmodified approach explained above.

    Let the modified feedback nonholonomic criteria be:

    (NH.1) iΞCHω=dH+λαCdΦα, for smooth functions λα:TQR, and

    (NH.2') LΞCHΦα=καΦα for some positive constants κα,

    where there is no summation in (NH.2').

    This approach has the advantage that if x0M, its trajectory exponentially approaches M, i.e. suppose that φt is the flow of the extended feedback nonholonomic vector field, then

    Φα(φt(x0))=eκαΦα(x0)0.

    Mimicking the derivation of (3.1), we arrive at

    iΞCHω=dHmαβ({H,Φα}καΦα)CdΦβ. (3.2)

    Example 1 (Simple Pendulum). Consider the simple pendulum viewed as a constrained system in R2. The Hamiltonian and constraint are

    H=12m(p2x+p2y)+mgy,Φ=xpx+ypy=0,

    where m is the mass of the pendulum and g is the acceleration due to gravity. The equations of motion from (3.2) are

    ˙x=1mpx,˙y=1mpy,˙px=xx2+y2[mgy1m(p2x+p2y)κ(xpx+ypy)],˙py=yx2+y2[mgy1m(p2x+p2y)κ(xpx+ypy)]mg. (3.3)

    A valid initial condition, z0=(x(0),y(0),px(0),py(0)), for (3.3) needs to satisfy the constraint Φ(z0)=0. However, due to running numerical errors, the constraint will not be preserved as the state is integrated. By choosing κ>0, a correcting feedback is introduced to help preserve Φ=0 as shown in Figure 1.

    Figure 1.  The introduction of a positive κ helps to stabilize the numerical trajectory. Numerical integration was performed with Matlab's ode45.

    We end this section with a brief consideration of ideal constraints - constraints that perform no work on the system. As discussed in [31], a constraint Ψ on TQ is ideal if and only if

    ˙qiΨ˙qi|M=0.

    On the cotangent side, this condition becomes

    CdΦ(XH)|M=0.

    In particular, if all the constraints are ideal, then energy is conserved. Indeed,

    ˙H=mαβ{H,Φα}CdΦβ(XH).

    Constraints will not be assumed to be ideal. Indeed, examples 7.1.2, 7.1.3, and 7.2.1 will all be subject to non-ideal constraints. In particular, these examples will be volume-preserving but not energy-preserving (although a modified energy integral exists for the first two examples [34]).

    For a given nonholonomic vector field XMHX(M), a volume-form μΩ2nk(M) is invariant if and only if

    divμ(XMH)=0.

    This section offers a way to choose a volume-form μ based off a (non-canonical choice) of realization C as well as showing how to compute the divergence of the nonholonomic vector fields with respect to this volume.

    The symplectic manifold TQ has a canonical volume form ωn. However, the nonholonomic flow takes place on a submanifold MTQ which is 2nk dimensional. Therefore, ωn is not a volume form on M. Here, we construct a volume form on M which is unique up to the choice of realization. The derivation of this will be similar to the construction of the volume form on an energy surface in §3.4 of [35]. For the realization C={Φ1,,Φk}, define the k-form

    σC:=dΦ1dΦk.

    Definition 4.1. If we denote the inclusion map by ι:MTQ, then a nonholonomic volume, μC, is given by

    μC=ιε,σCε=ωn.

    Proposition 4.2. Given an ordered collection of constraints, C, the induced volume form μC is unique.

    Proof. Suppose that ε and ε are two forms satisfying σCε=ωn. Then

    εε=α,σCα=0.

    Now let ι:MTQ be the inclusion. Then from the above, we see that

    ιε=ιε+ια.

    The result will follow so long as ια=0. Suppose that ια0 and choose vectors v1,,v2nkTxMTxTQ such that α(v1,,v2nk)0. Complete this collection of vectors to a basis of TxTQ: v1,,v2nk,v2nk+1,,v2n such that σC(v2nk+1,,v2n)0. Then we have

    σCα(v1,,v2n)=(1)(2nk)kα(v1,,v2nk)σC(v2nk+1,,v2n)0,

    which is a contradiction.

    Remark 4.3. Notice that for an ordered collection of constraints the volume form is unique. However, changing the order of the constraints changes the sign of the induced volume form and re-scaling constraints re-scales the volume form. In this sense, C uniquely determines μC, but M only determines μC up to a multiple. This is to be expected as all volume-forms are related by a multiplicative factor.

    While examining the failure of Liouville's theorem (Theorem 2.2) for nonholonomic systems, we will see when μC is preserved under the flow of XMH. More generally, we will consider the existence of a smooth density fC(M) when fμC is preserved.

    Let ω=dqidpi be the standard symplectic form on TQ. This in turn induces a volume form ωn. A measure of how much a flow fails to preserve a volume form is described by its divergence. We proceed with computing the divergence of a nonholonomic vector field, divμC(XMH). When this is nonzero, we will be interested in finding a density, f, such that divfμC(XMH)=0. This problem will be addressed in §5.

    Before we begin with the divergence calculation, we first present a helpful lemma which allows us to relate the divergence of the extended nonholonomic vector field with the corresponding restricted vector field.

    Lemma 4.4. If C is a realization of the constraint manifold MTQ, then

    divωn(ΞCH)|M=divμC(XMH).

    Proof. Leibniz's rule for the Lie derivative provides

    LΞCHωn=LΞCH(σCε)=(LΞCHσC)ε+σC(LΞCHε).

    However, LΞCHσC=0 because the constraints are preserved under the flow. Applying this, we see that

    LΞCHωn=σC(LΞCHε),

    which gives

    (divωn(ΞCH))σCε=σC(LΞCHε).

    Due to the fact that the Lie derivative commutes with restriction, the result follows.

    This lemma allows for us to calculate the divergence of the extended nonholonomic vector field and to restrict to the constraint distribution afterwards.

    Using Cartan's magic formula along with Lemma 4.4, we see

    LΞCH(ωn)=iΞCHdωn+diΞCHωn=nd(iΞCHωωn1)=nd(iΞCHω)ωn1n(iΞCHω)dωn1=ndνCHωn1.

    The divergence of the system is controlled by the failure of νCH to be closed. Its derivative is given by

    dνCH=ddHd(mαβ{H,Φα}CdΦβ).

    Let ϕα=mαβΦβ. As Φβ=0 on M, the derivative becomes (where restriction to M is implied)

    dνCH=d{H,ϕβ}CdΦβ+{H,ϕβ}dCdΦβ.

    The form dνCH is a general 2-form. However, the only terms that survive once it is wedged with ωn1 are the diagonal terms, dqidpi. Using local coordinates, aβidqi=CdΦβ, we have

    (dνCH)diag=({H,ϕβ}piaβi{H,ϕβ}aβipi)dqidpi.

    Therefore, the divergence is given by

    divμC(XMH)=n({H,ϕβ}piaβi{H,ϕβ}aβipi).

    The first component of this expression can be written intrinsically. Notice that

    dπQXf=fpiqi,

    where πQ:TQQ is the cotangent bundle projection. Therefore, the divergence can now be written as

    divμC(XMH)=n(CdΦβ([XH,Xϕβ])+{H,ϕβ}aβipi). (4.1)

    Suppose the Lagrangian is natural with Riemannian metric g=(gij). Then

    CdΦβ=gijΦβpjdxi,aβi=gijΦβpj,aβipi=gij2Φβpipj.

    Let us call this final double sum Mβ. Applying this to (4.1), we have

    divμC(XMH)=nCdΦβ([XH,Xϕβ])n{H,ϕβ}Mβ.

    Next, assume that (in addition to the Lagrangian being natural), the constraints are affine which makes Mβ=0. Also the matrix, mαβ, does not depend on p and therefore, the divergence can be written as

    divμC(XMH)=nmαβπQηβ([XH,XΦα]), (4.2)

    as CdΦβ=πQηβ for affine constraints.

    For affine systems with natural Lagrangians, the divergence can be expressed intrinsically via (4.2). To understand the divergence in the general nonlinear case (with potentially non-natural Lagrangians), we need an intrinsic way to interpret the term Mβ. This can be accomplished via the following procedure:

    The divergence for a nonholonomic system subject to nonlinear constraints is thus

    divμC(XMH)=dim(Q)(CdΦβ([XH,Xϕβ])+{H,ϕβ}ΔC(Φβ)), (4.3)

    which is the key to proving Theorem 5.2 below.

    In general, the divergence of a nonholonomic system does not vanish as (4.1) shows. When does there exist a different volume form on M that is invariant under the flow? i.e. does there exist a density ρ>0 such that divρμC(XMH)=0? Finding such a ρ requires solving a certain type of partial differential equation which is known as the smooth dynamical cohomology equation. Solving this PDE is generally quite difficult. However, when the constraints are affine and the density is assumed to be basic, the problem becomes considerably more tractable

    What conditions need to be met for ρ such that ρμC is an invariant volume form? Using the formula for the divergence as well as the fact that the Lie derivative is a derivation yields:

    divρμC(XMH)=divμC(XMH)+1ρLXMH(ρ).

    Therefore the density, ρ, yields an invariant measure if and only if

    1ρLXMH(ρ)=divμC(XMH). (5.1)

    Notice that the left hand side of (5.1) can be integrated to

    1ρLXMH(ρ)=d(lnρ)(XMH).

    Calling λ=lnρ, we have the following lemma.

    Lemma 5.1. For a nonholonomic vector field, XMH, there exists a smooth invariant volume, ρμC, if there exists an exact 1-form α=dλ such that

    α(XMH)=divμC(XMH). (5.2)

    Then the density is (up to a multiplicative constant) ρ=eλ.

    Therefore the existence of invariant volumes boils down to finding global solutions to the PDE (5.2). Using the divergence calculation, (4.3), we can state the main result.

    Theorem 5.2. The nonholonomic Hamiltonian vector field XMH possesses a smooth invariant volume, ρμC, if there exists an exact 1-form αΩ1(TQ) such that

    α(XMH)=dim(Q)(CdΦβ([XH,Xϕβ])+{H,ϕβ}ΔC(Φβ)).

    The remainder of this section deals with uniqueness of solutions and a necessary condition for solutions to exist.

    Remark 5.3. PDEs of the form dg(X)=f for a given smooth function f, vector field X and with g as the unknown are called cohomology equations [36,37]. Thus the equation (5.2) is a cohomology equation.

    The problem of existence is quite difficult in general and we postpone that discussion until the next subsection where we assume that the solution has the form ρ=πQ˜ρ. In the meantime, assuming that there exists a function λC(M) that solves (5.2), do there exist other solutions? Suppose that λ1 and λ2 both solve (5.2). Then their difference must be a first integral of the system: LXMH(λ1λ2)=0. Solutions of (5.2) are then unique up to constants of motion. i.e. if λ solves (5.2), then every invariant density has the form (again, up to a multiplicative constant)

    ρ=exp(λ+constantofmotion).

    Therefore invariant measures can be thought of as an affine space with dimension being equal to the number of first integrals of the nonholonomic system.

    In general, solving the cohomology equation (5.2) is quite difficult. It turns out, however, that it is relatively easy to determine necessary and sufficient conditions on the solvability when the density is assumed to be basic and the constraints are affine.

    For the remainder of this section, unless otherwise stated, the constraints are assumed to be affine. In particular, if NTQ, then the constraints have the form

    Ψα(v)=ηα(v)+τQξα(v),

    for 1-forms ηα and function ξα. On the cotangent side, the constraints become

    Φα(p)=P(Wα)(p)+πQξα(p), (5.3)

    where Wα=FL1ηα.

    Definition 5.4. A density ρ:TQR is said to be basic if ρ=πQ˜ρ for some ˜ρ:QR.

    Under this assumption, (4.2) can be presented in a surprisingly nice way. In this case, the divergence can be described by an equivalence class of "affine-forms." The density form and density function, defined below, is a representative element from this class.

    Definition 5.5. Let C be an affine realization of MTQ of the form (5.3). Then, define the density form and the density function to be

    ϑC=mαβLWαηβ,ζC=mαβLWαξβ,

    respectively.

    Studying this pair, (ϑC,ζC), provides necessary and sufficient conditions for the existence of basic densities.

    Theorem 5.6. For a natural Hamiltonian system subject to affine constraints, there exists an invariant volume of the form (πQ˜ρ)μC if and only if there exists functions φγ such that

    nϑCφγηγ=dln˜ρ,nζCφγξγ=0. (5.4)

    Proof. To prove this, we will show that nπQϑ(XMH)nπQζC=divC(XMH). Recall that the differential of a 1-form is given by dα(X,Y)=Xα(Y)Yα(X)α([X,Y]) and that πQηβ(XH)|M=ξβ. Returning to (4.2), we have

    divμC(XMH)=nmαβπQηβ([XH,XΦα])=nmαβ(XHmαβ+XΦαξβπQdηα(XH,XΦα))=nmαβ(dmαβ(˙q)+dξβ(Wα)dηα(˙q,Wβ))=nmαβ[(diWβηα+iWβdηα)(˙q)+dξβ(Wα)]=nmαβ[LWβηα(˙q)+LWαξβ]=nϑC(˙q)nζC.

    This computation shows that divμC(XMH)=nϑC(˙q)nζC, but ˙q cannot be arbitrary as it must lie within M which states that ηα(˙q)+ξα=0. Adding multiples of the constraints to the divergence to produce an exact 1-form yields (5.4).

    This theorem allows for a straight-forward algorithm to find invariant volumes in nonholonomic systems subject to affine constraints; one only needs to compute the pair (ϑC,ζC) and determine whether or not it can be made exact by appending constraints to it. This procedure will be carried out on multiple examples in §7.

    Remark 5.7. In the pure kinetic energy case with linear constraints discussed in [24], it is proved that if the system admits an (arbitrary) invariant volume, then one can always find another invariant volume form whose density function depends only on the (reduced) configuration variables. Moreover, systems subjected to affine constraints may possess an invariant volume whose density is not basic, cf. [38].

    The above shows that "exactness" of (ϑC,ζC) determines the existence of a density depending on configuration. How does this depend on the choice of C to realize the constraints? It turns out the answer is independent of the choice of realization.

    Theorem 5.8. Let C1 and C2 both be affine realizations of the constraint manifold M. Suppose that there exist functions φγ such that

    ϑC1φγηγ1=dln˜ρ1,ζC1φγξγ1=0,

    then there exists functions ψγ such that

    ϑC2ψγηγ2=dln˜ρ2,ζC2ψγξγ2=0.

    Moreover, (πQ˜ρ1)μC1=(πQ˜ρ2)μC2 modulo a constant of motion.

    Proof. Let ηγ2=cγαηα1 and ξγ2=cγαξα1 where cγα is the required coordinate change. For ease of notation, we will drop the subscript "1." The matrices are related via

    mαβ2=cαγcβδmγδ1.

    The density form is

    mαβ2ϑ2=LcαγWγ(cβδηδ)=cαγcβδLWγηδ+cβδmδγ1dcαγ+cαγdcβδ(Wγ)ηδ=mαβ2ϑ1+mαβ2σγαdcαγ+cαγdcβδ(Wγ)ηδ,

    where (σγα):=(cαγ)1. By Lemma 5.9 below, we have

    mαβ2ϑ2=mαβ2ϑ1+mαβ2d[lndet(cαγ)]+cαγdcβδ(Wγ)ηδ. (5.5)

    This shows that ϑ1 and ϑ2 differ by something exact and a multiple of the constraining 1-forms. Using the fact that ϑ1φγηγ1=dln˜ρ1, we see that

    ϑ2(φγ+Cγ)ηγ=d[ln˜ρ1+lndet(cαγ)],

    where mαβ2Cγ=cαδdcβγ(Wδ). This provides

    ˜ρ2=det(σγα)˜ρ1,ψγ=cδγ[φδ+Cδ].

    These values of ψγ make ϑ2 exact. We next check that these values of ψγ make ζ2 vanish.

    mαβ2(ζ2ψγcγδξδ)=cαγcβδLWγξδ+cαγdcβδ(Wγ)ξδmαβ2ψγcγδξδ=mαβ2ζ1mαβ2φδξδ=0.

    It remains to show that (πQ˜ρ1)μC1=(πQ˜ρ2)μC2. This is equivalent to μC1=(πQdet(cαγ))μC2. Recalling Definition 4.1, we have σ1=dΦ1dΦk and

    σ2=ki=1d(ciγΦγ)=det(cαγ)σ1+ki=1Φγdciγ.

    Upon applying the constraints, the last term disappears and we obtain

    σ1ε1=σ2ε2=det(cαγ)σ1ε2.

    This implies the desired result.

    A reason why studying (ϑC,ζC) is insightful is that it immediately demonstrates why holonomic systems systems are measure-preserving, i.e. when ηα=dfα and ξα=0, (5.4) can always be solved. This short section idemonstrates how our general theory collapses to the known holonomic case. We start with a helpful lemma.

    Lemma 5.9 (*). mαβdmαβ=d[lndet(mαβ)].

    *We thank Dr. Alexander Barvinok for help with this proof.

    Proof. It suffices to check along a curve in the manifold. Let γ:IQ be a curve and let A(t)=(mαβ)γ(t) be the mass matrix along the curve. Note that A(t) is positive-definite and changes smoothly with t. We have

    ddtlndetA(t)=ddtdetA(t)detA(t)=mi=1detAi(t)detA(t),

    where Ai(t) is obtained from A(t) by differentiating the i-th row and leaving all other rows intact, i.e.

    Ai(t)=(a11(t)a1m(t)a(i1)1(t)a(i1)m(t)ai1(t)aim(t)a(i+1)1(t)a(i+1)m(t)am1(t)amm(t)).

    Expanding detAi(t) along the i-th row:

    detAi(t)=mj=1(1)i+j1aij(t)detAij(t),

    where Aij(t) is the (m1)×(m1) matrix obtained from Ai(t) and hence from A(t) by crossing out the i-th row and j-th column.

    Next, observe that (1)i+j1detAij/detA(t) is the (j,i)-th entry of the inverse matrix A1(t)=(bij)(t), and since A(t) is symmetric, is also the (i,j)-th entry of (bij)(t). Summarizing,

    ddtlndetA(t)=mi,j=1aij(t)bij(t).

    Proposition 5.10. If the constraints are holonomic, then there exist function φγ such that ϑCφγηγ is exact. In particular, if C is chosen such that all ηα are closed, ϑC is exact.

    Proof. When the constraints are holonomic, the 1-forms ηα can be chosen such that they are closed. Then the density form is

    which is exact by Lemma 5.9. If a different realization is chosen, Theorem 5.8 states that the resulting (5.4) is still solvable.

    It turns out that for natural Lagrangian systems subject to linear constraints, the divergence - particularly the density form - is encoded in the nonholonomic connection. This interpretation seems to be new.

    Throughout this section, let L:TQR be a natural Lagrangian with Riemannian metric g subject to the linear constraints ηα(v)=0. The nonholonomic connection for this system is given by (cf. §5.3 in [20] and [39]):

    CXY=XY+Wimij[X(ηj(Y))ηj(XY)].

    The equations of motion can then be described via

    C˙q˙q=F,

    where F contains the potential and external forces (the constraint forces are contained in the connection).

    The nonintegrability of the constraints appears in the torsion of the connection. Computing this, we see

    TC(X,Y)=CXYCYX[X,Y]=Wimij[X(ηj(Y))Y(ηj(X))ηj(XYYX)]=Wjmij[X(ηj(Y))Y(ηj(X))ηj([X,Y])]=Wjmijdηj(X,Y).

    The torsion can be written as

    TC=mαβWαdηβ.

    Indeed, if the constraining 1-forms ηj are all closed (so holonomic) then the torsion vanishes. It is worth pointing out that the torsion is vertical-valued; if X,YM, then TC(X,Y)M, i.e. T(X,Y) is orthogonal to the constraint distribution.

    Due to the fact that the torsion is a (1, 2)-tensor, its trace will be a (0, 1)-tensor. Therefore, the trace of the nonholonomic torsion will be a 1-form:

    trTC=mαβiWαdηβ.

    Returning to the density form, we see that

    trTC+dlndet(mαβ)=ϑC,

    i.e. the trace of the torsion differs from the density form by something exact. This leads to the following theorem.

    Theorem 6.1. A natural nonholonomic system subject to linear constraints has an invariant volume of the form (πQρ)μC if and only if there exists functions φγ such that

    trTC+φγηγ

    is exact.

    Remark 6.2. The vanishing of the torsion shows that the constraints are integrable while the integrability of the (trace of the) torsion shows that a volume is preserved.

    In the case of nonholonomic systems, the nonholonomic connection is compatible with the metric but has nonzero torsion. This idea extends to arbitrary, metric-compatible connections as the following theorem states.

    Theorem 6.3. Let ˜ be an affine connection compatible with the metric with torsion ˜T. There exists an invariant volume with density of the form πQρ for the geodesic spray if and only if tr˜T is exact.

    Proof. Consider the volume form on TQ given by

    Ω=detgdx1dxndv1dvn.

    We want to compute LXΩ where X is the geodesic spray given by

    X=vixiΓijkvjvkvi.

    The Lie derivative is then

    LXΩ=diXΩ=(d[detg](v)detg(Γiik+Γiki)vk)1detgΩ,

    and therefore the divergence is given by

    divΩ(X)=d[lndetg](v)(Γiik+Γiki)vk. (6.1)

    We will now use the fact that the connection is compatible with the metric:

    gjkxi=gkΓij+gjΓik.

    This implies that

    gjkgjkxi=gjkgkΓij+gjkgjΓik=δjΓij+δkΓij=2Γkik.

    Integrating the left-hand side above gives

    d[lndetg](v)=2Γikivk. (6.2)

    Substituting (6.2) into (6.1), we get

    divΩ(X)=(ΓikiΓiik)vk.

    It remains to show that this is the trace of the torsion. Indeed,

    ˜T=(ΓkijΓkji)xkdxidxjtr˜T=(ΓikiΓiik)dxk.

    We conclude that

    divΩ(X)=tr˜T(v).

    This shows that a way to interpret the torsion of a connection is by measuring how much the geodesic spray fails to preserve volume.

    We end this work with applying both Theorem 5.2 and Theorem 5.6 to various nonholonomic systems by calculating the general divergence for nonlinear systems or ϑC and ζC for affine/linear systems. Examples are taken from [18,20,31].

    We begin by examining multiple nonholonomic systems subject to affine/linear constraints.

    As an example of Theorem 5.6, we will prove that no invariant basic volumes exist for the Chaplygin sleigh. The Chaplygin sleigh has the configuration space Q=SE2, the special Euclidean group, and the following Lagrangian

    L=12(m˙x2+m˙y2+(I+ma2)˙θ22ma˙x˙θsinθ+2ma˙y˙θcosθ),

    where (x,y)R2 is the coordinate of the contact point, θSO2 is its orientation, m is the sleigh's mass, I is the moment of inertia about the center of mass, and a is the distance from the center of mass to the contact point (cf. §1.7 in [20]).

    The nonholonomic constraint is that the sleigh can only slide in the direction it is pointing and is given by

    ˙ycosθ˙xsinθ=0,

    which corresponds to the 1-form η=(cosθ)dy(sinθ)dx and function ξ=0.

    To determine the existence of an invariant volume, we only need to compute ϑC as ζC=0. The constraining vector field and 1-form are:

    W=ma2+IIm[cosθysinθx]aIθ,η=(cosθ)dy(sinθ)dx.

    This gives us

    ϑC=1η(W)LWη=mama2+I[(sinθ)dy+(cosθ)dx].

    As a consequence of this, the divergence of the Chaplygin sleigh is given by

    divμC(XDH)=3mavI+ma2,v=˙xcosθ+˙ysinθ. (7.1)

    We want to show that for any ˜ηΓ(D0), ϑC+˜η is not exact. Because there is only one constraint, it suffices to show that there does not exist a smooth k such that ϑC+kη is closed, i.e. it requires the following to be zero:

    d(ϑC+kη)=mama2+I[(cosθ)dθdy(sinθ)dθdx]+(kxcosθ+kysinθ)dxdy+(kθcosθksinθ)dθdy(kθsinθ+kcosθ)dθdx.

    Separating the above, we need the following three to vanish:

    0=kxcosθ+kysinθ,0=kθcosθksinθ+mama2+Icosθ,0=kθsinθ+kcosθ+mama2+Isinθ. (7.2)

    The second two lines of (7.2) are overdetermined for k in the θ-direction and are inconsistent (unless a=0 and we obtain the trivial solution k0). Therefore, there does not exist a smooth k such that ϑC+kη is closed. We note that this is compatable with the known result that when a=0, no asymptotically stable dynamics occur.

    The next example is that of the falling rolling disk whose configuration space is Q=SE2×S1. Its Lagrangian is

    L=m2[(ξR(˙φsinθ+˙ψ))2+η2sin2θ+(ηcosθ+R˙θ)2],

    where

    ξ=˙xcosφ+˙ysinφ+R˙ψ,η=˙xsinφ+˙ycosφ.

    We will consider both the case of the rolling coin on a stationary table (linear constraints) and on a rotating table (affine constraints).

    Stationary Table The case of the disk on a stationary table was studied in [40]. The constraints are given by the vanishing of the following 1-forms:

    η1=cosφdx+sinφdy+Rdψ,η2=sinφdx+cosφdy. (7.3)

    The corresponding dual vector fields are

    W1=1mcosφx+1msinφy+RIψ,W2=J+mR2Jm+m2R2sin2θ[sinφx+cosφy]RcosθJm+m2R2sin2θθ.

    Computing ϑC, we obtain

    ϑC=mR2sin(2θ)J+mR2sin2θ,

    which is exact. The resulting invariant volume is

    1J+mR2sin2θμC.

    Rotating Table Although the case of the stationary table has been studied, the authors are unaware of any results for the case of a rotating table. If the falling rolling disk is placed on a table with constant angular velocity Ω, the constraints become affine with

    η1(v)+ξ1=0,ξ1=Ω(ycosφxsinφ)η2(v)+ξ2=0,ξ2=Ω(xcosφ+ysinφ)

    where ηα are from (7.3). The volume from the stationary table is still preserved as

    LW1ξ1=LW2ξ2=0.

    We next consider the case of a non-homogeneous sphere rolling without slipping on a horizontal plane, both stationary and rotating. The center of mass of the sphere is located at its geometric center while its principal moments of inertia are distinct. This example has been studied by Chaplygin [41] and, e.g. [1,42,43].

    The Lagrangian is the kinetic energy,

    L=12I1(˙θcosψ+˙φsinψsinθ)2+12I2(˙θsinψ+˙φcosψsinθ)2+12I3(˙ψ+˙φcosθ)2+12M(˙x2+˙y2),

    where it is assumed that the radius is 1.

    Stationary Table When the table is stationary, the constraints are given by the two 1-forms

    η1=dxsinφdθ+cosφsinθdψ,η2=dy+cosφdθ+sinφsinθdψ.

    The density form is given by ϑC=Adθ+Bdψ where

    A=Msin(2θ)[J1+J2sin2ψ]2(J3+J4sin2θ+J5sin2θsin2ψ),B=MJ2sin(2ψ)sin2θ2(J3+J4sin2θ+J5sin2θsin2ψ),

    and the constants Jj are

    J1=I1I2I1I3+I2MI3M,J2=I1I3I2I3+I1MI2M,J3=I3M2+I1I2I3+I1I3M+I2I3M,J4=I2M2I3M2+I1I2MI1I3M,J5=I1M2I2M2+I1I3MI2I3M.

    The density form is exact and the resulting invariant volume is

    (1+β1β)12cosθ(J3+(J3+J4)tan2θ)12μC,β=J5sin2θsin2ψ2J3+2J4sin2θ+J5sin2θsin2ψ.

    Rotating Table When the table is rotating, the constraints becomes affine

    η1(v)+ξ1=0,ξ1=Ωy,η2(v)+ξ2=0,ξ2=Ωx.

    Notice that,

    LW1ξ1=LW2ξ2=0,LW2ξ1=LW1ξ2=ΩM.

    As the matrix (mαβ) is symmetric, we have that the product

    mαβLWαξβ=0.

    Therefore, the Chaplygin sphere on a rotating table is volume-preserving with the same volume as in the stationary case.

    Theorem 5.2 does not require that Q be orientable. Consider the Möbius strip immersed in R3 by

    x=(1+vcos(u2))cos(u),y=(1+vcos(u2))sin(u),z=vsin(u2), (7.4)

    for 0u2π and 1/2<v<1/2. The Euclidean metric pulled back to the Möbius strip is

    g=(4vcos(u2)+2v2cos(u2)+v24+2)dudu+2dvdv.

    Notice that the above metric is for the double cover rather than the Möbius strip itself as a consequence of (7.4) only being an immersion. However, the resulting equations of motion will be well-defined on the strip. In two-dimensional systems, any single constraint is automatically holonomic which will make volume-preservation trivial. Let us "thicken" the strip by w with resulting metric

    gthick=(4vcos(u2)+2v2cos(u2)+v24+2)dudu+2dvdv+dwdw.

    Consider the linear nonholonomic constraint

    η=dv+sin(u)dw,W=12v+12sin(u)w.

    The density form is

    ϑC=sin(u)cos(u)1+sin2(u)du,

    which is exact. This produces the invariant volume

    1+sin2(u)μC,

    which is shown in Figure 2.

    Figure 2.  A plot of the density corresponding to the invariant volume for the nonholonomic system on the Möbius strip.

    The previous examples consisted of a natural Lagrangian and affine constraints which were describable by Theorem 5.6. The next two examples utilize Theorem 5.2 on systems with natural Lagrangians but nonlinear constraints. These examples can be found in [31,44].

    Let L:TR3R be the Lagrangian given by

    L=12m(˙x2+˙y2+˙z2)mgz,

    subject to the nonlinear constraint of constant kinetic energy

    Ψ=˙x2+˙y2+˙z2c=0,c>0.

    Transferring to the Hamiltonian side, we have

    H=12m(p2x+p2y+p2z)+mgz,Φ=12(p2x+p2y+p2z)c=0,

    where Φ is normalized. As the constraints are nonlinear, (4.3) can be used to determine the divergence. The requisite data to compute the divergence is

    CdΦ=m(pxdx+pydy+pzdy),CdΦ(XΦ)=m(p2x+p2y+p2z),ΔCΦ=3m,ϕ=12mcm(p2x+p2y+p2z),{H,ϕ}=2gcpz(p2x+p2y+p2z).

    Although CdΦ(XΦ)0 everywhere, it does in some tubular neighborhood of M. The divergence of the system is given by

    divμC(XMH)=3CdΦ([XH,Xϕ])3{H,ϕ}M=18mgcpz(p2x+p2y+p2z)218mgcpz(p2x+p2y+p2z)2=9mgpzc,

    where we used the fact that p2x+p2y+p2z=2c. The divergence does not vanish so μC is not preserved. However, there exists an exact 1-form which produces the divergence:

    9mgcdz(XMH)=divμC(XMH).

    Therefore, the following volume-form is preserved

    exp(9mgzc)μC.

    The other nonlinear constraint example we will examine is Appel's example. The Lagrangian is the same as the constant kinetic energy case,

    L=12m(˙x2+˙y2+˙z2)mgz,

    while the nonlinear constraint is now

    Ψ=a2(˙x2+˙y2)˙z2=0.

    The data on the Hamiltonian side becomes

    H=12m(p2x+p2y+p2z)+mgz,Φ=a22(p2x+p2y)12p2z=0,

    where, again, the constraint is normalized. As before, the constraints are nonlinear so (4.3) can be used. The requisite data is

    CdΦ=a2mpxdx+a2mpydympzdz,CdΦ(XΦ)=ma4(p2x+p2y)+mp2z,ΔCΦ=m(2a21),ϕ=a2(p2x+p2y)p2z2ma4(p2x+p2y)+2mp2z,{H,ϕ}=a2(a2+1)gpz(p2x+p2y)(a4p2x+a4p2y+p2z)2.

    The constraint is not admissible as CdΦ(XΦ) vanishes at px=py=pz=0 which is in the constraint manifold. As long as a1, this is the only place where this degeneracy occurs. The divergence is

    divμC(XMH)=3CdΦ([XH,Xϕ])3{H,ϕ}M=12a2(a2+1)mgpz(p2x+p2y)(a6p2x+a6p2yp2z)(a4p2x+a4p2y+p2z)3.

    To simplify the divergence, notice that the constraint makes p2z=a2(p2x+p2y). Substituting this, the divergence becomes

    divμC(XMH)=12(a21)mg(a2+1)pz.

    The following exact 1-form solves the cohomology equation

    12(a21)a2pzdpz(XMH)=divμC(XMH),

    as

    ˙pz=a2mg1+a2.

    Therefore, the following form is preserved

    pKzμC,K=12(a21)a2.

    Unfortunately, this form is not a volume-form as it vanishes when pz=0.

    All of the examples examined so far have been for systems whose Lagrangian is natural. We conclude this work with two examples whose Lagrangian is not given by a Riemannian metric.

    Consider the Lagrangian

    L=14(˙x4+˙y4+˙z4),

    subject to the nonintegrable constraint ˙z=x˙y. The Hamiltonian is

    H=34(p4/3x+p4/3y+p4/3z),

    and the constraint becomes

    Φ=p1/3zxp1/3y.

    Notice that singularities appear when the momentum vanishes. Continuing with the computations, we have

    CdΦ=dzxdyCdΦ(XΦ)=x23p2/3y+13p2/3zϕ=3pz3x2pz+3xpyx2+ν2/3,ν=pypzΔCΦ=0

    Applying the constraint, the divergence simplifies to

    divμC(XMH)=3p1/3xx42x(1+x4)=3x42x(1+x4)dx(XMH)

    The following form is exact

    α=3x46x(1+x4)dx,

    and an invariant form is

    (1+x4)9/4x6μC.

    Notice the singularity at x=0.

    Suppose we have the relativistic Lagrangian

    L=m0c21v2c2,v2=˙x2+˙y2+˙z2.

    With the same constraints ˙z=x˙y. The Hamiltonian and constraint are

    H=cp2+m2c2,Φ=c(pzxpy)p2+m2c2.

    The computation yields

    CdΦ=dzxdyCdΦ(XΦ)=c((m2c2+p2x)(1+x2)+(xpz+py)2)(m2c2+p2)3/2ϕ=(p2+m2c2)(pzxpy)(m2c2+p2x)(1+x2)+(xpz+py)2ΔCΦ=0

    Applying the constraints, the divergence simplifies to

    divμC(XMH)=3xpx(1+x2)m2c2+p2x+(1+x2)p2y=3x1+x2dx(XMH)

    Therefore, an invariant volume is given by

    (1+x2)3/2μC.

    We thank Dr. J.C. Marrero for pointing us to existing work in this field as well as the reviewers for helpful comments.



    [1] Z. W. Shen, H. Z. Zhou, H. Yin, J. T. Wu, L. Li, Fine adjusted-customized ablation LASIK treatment for myopia, Int. J. Ophthalmol., 5 (2005), 1194-1197.
    [2] Z. J. Fan, S. J. Xu, Z. H. Jia, B. C. Liu, Clinical significance of corneal Q value in myopic patients, Int. J. Ophthalmol., 6 (2006), 642-643.
    [3] X. U. Li, Q. Tao, L. Y. Zhuang, Visual outcome after optimized aspheric transitionzone laser keratomileusis compared toconventional LASIK, Int. Eye Sci., 7 (2007), 623-625.
    [4] S. MacRae, Excimer ablation design and elliptical transition zones, J. Cataract Refractive Surg., 25 (1999), 1191-1197. doi: 10.1016/S0886-3350(99)00144-3
    [5] L. Fang, Y. Wang, X. He, Theoretical analysis of wavefront aberration caused by treatment decentration and transition zone after custom myopic laser refractive surgery, J. Cataract Refractive Surg., 39 (2013), 1336-1347. doi: 10.1016/j.jcrs.2013.03.020
    [6] G. M. Dai, Application of complementary error function to the transition zone of myopic ablation shapes in refractive surgery, J. Appl. Math. Phys., 5 (2017), 1521-1528. doi: 10.4236/jamp.2017.58125
    [7] E. Gross, R. Hofer, J. Wong, Application of blend zones, depth reduction, and transition zones to ablation shapes, Patent No. 8216213, 2012.
    [8] J. Shen, Y. Zhang, W. Liao, Mathematical model based corneal toric surface for excimer laser refractive surgery, J. Southeast Univ., 36 (2006), 531-536.
    [9] M. J. Endl, C. E. Martinez, S. D. Klyce, M. B. McDonald, S. J. Coorpender, R. A. Applegate, et al., Effect of larger ablation zone and transition zone on corneal optical aberrations after photorefractive keratectomy, Arch. Ophthalmol., 119 (2001), 1159-1164. doi: 10.1001/archopht.119.8.1159
    [10] I. B. Damgaard, M. Ang, A. M. Mahmoud, M. Farook, C. J. Roberts, J. S. Mehta, Functional optical zone and centration following SMILE and LASIK: a prospective, randomized, contralateral eye study, J. Refractive Surg., 35 (2019), 230-237.
    [11] A. P. Nisarta, D. Desai, K. Solanki, Corneal higher order aberrations after aspheric LASIK treatment, Int. J. Health Sci. Res., 6 (2016), 107-113.
    [12] T. Gamaly, LASIK with the optimized aspheric transition zone and cross-cylinder technique for the treatment of astigmatism from 1.00 to 4.25 diopters, J. Refractive Surg., 25 (2009), S927-930.
    [13] Y. H. Komai, I. Toda, N. A. Kato, M. Ito, T. Yamaoto, K. Tsubota, Comparison of LASIK using the NIDEK EC-5000 optimized aspheric transition zone (OATz) and conventional ablation profile, J. Refractive Surg., 22 (2006), 546-555. doi: 10.3928/1081-597X-20060601-06
    [14] M. Hantera, Comparison of postoperative wavefront aberrations after NIDEK CXIII optimized aspheric transition zone treatment and OPD-guided custom aspheric treatment, J. Refractive Surg., 25 (2009), S922-926.
    [15] M. Balidis, Biomechanical profile of refractive surgery procedures, Acta Ophthalmol., 97 (2019).
    [16] G. V. Voronin, I. A. Bubnova, Changes in biomechanical properties of the cornea after keratorefractive surgery, Vestn. Oftalmologii, 135 (2019), 108-112.
    [17] D. V. Franus, Change in the stress-strain state of the cornea after refractive surgery, in International Conference on Mechanics-seventh Polyakhovs Reading, IEEE, (2015), 1-4.
    [18] M. A. Widlicka, W. Srodka, P. K. Berkowska, The biomechanical modelling of the refractive surgery, Optik, 120 (2009), 923-933. doi: 10.1016/j.ijleo.2008.03.026
    [19] I. Simonini, A. Pandolfi, Customized finite element modelling of the human cornea, PloS one, 10 (2015), e0130426.
    [20] K. Salmenhaara, Impact of Refractive Surgery to the Biomechanical Properties of the Cornea-a Finite Element Analysis, Master thesis, Aalto University, 2015.
    [21] H. Y. Tian, B. Wang, J. X. Huo, Research progress of the effects of high order aberrations on visual quality after the LASIK, Prog. Mod. Biomed., 14 (2014), 3393-3395.
    [22] P. Vinciguerra, F. I. Camesasca, Treatment of hyperopia: a new ablation profile to reduce corneal eccentricity, J. Refractive Surg., 18 (2002), 315-317.
    [23] D. Epstein, P. Vinciguerra, M. Azzolini, P. Radiee, Long-term follow-up of hyperopic photorefractive keratectomy (prk) performed with a new algorithm, Invest. Ophthalmol. Visual Sci., 38 (1997).
    [24] Z. X. Dong, X. T. Zhou, Advences in biomechanical effects of laser corneal refractive surgery, Chin. J. Opthalmol., 48 (2012), 1053-1056.
    [25] J. Kim, Analysis of corneal higher-order aberrations after myopic refractive surgery, Curr. Opt. Photonics, 3 (2019), 72-77.
    [26] E. Lanchares, B. Calvo, M. A. D. Buey, J. A. Cristóbal, M. doblaré, The effect of intraocular pressure on the outcome of myopic photorefractive keratectomy: a numerical approach, J. Healthcare Eng., 1 (2010), 461-476.
    [27] L. Fang, W. Ma, Y. Wang, Y. Dai, Z. Fang, Theoretical analysis of wave-front aberrations induced from conventional laser refractive surgery in a biomechanical finite element model, Invest. Ophthalmol. Visual Sci., 61 (2020), 33-34.
  • This article has been cited by:

    1. Luis C. García-Naranjo, Rafael Ortega, Antonio J. Ureña, Invariant Measures as Obstructions to Attractors in Dynamical Systems and Their Role in Nonholonomic Mechanics, 2024, 29, 1560-3547, 751, 10.1134/S156035472456003X
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3092) PDF downloads(111) Cited by(1)

Figures and Tables

Figures(6)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog