# A symmetry and Noether charge preserving discretization of initial value
problems

ARAlexander Rothkopf  JNJan Nordström  Faculty of Science and Technology,
University of Stavanger, 4021, Stavanger, Norway  Department of Mathematics,
Linköping University, SE-581 83, Linköping, Sweden  Department of Mathematics
and Applied Mathematics, University of Johannesburg, P.O. Box 524, Auckland
Park 2006, Johannesburg, South Africa

###### Abstract

Taking insight from the theory of general relativity, where space and time are
treated on the same footing, we develop a novel geometric variational
discretization for second order initial value problems (IVPs). By discretizing
the dynamics along a world-line parameter, instead of physical time directly,
we retain manifest translation symmetry and conservation of the associated
continuum Noether charge. A non-equidistant time discretization emerges
dynamically, realizing a form of automatic adaptive mesh refinement (AMR),
guided by the system symmetries. Using appropriately regularized summation by
parts finite difference operators, the continuum Noether charge, defined via
the Killing vector associated with translation symmetry, is shown to be
exactly preserved in the interior of the simulated time interval. The
convergence properties of the approach are demonstrated with two explicit
examples.

Initial Value Problem, Summation By Parts, Time-Translation Invariance,
Conserved Noether Charge, Adaptive Mesh Refinement,

###### keywords:

PreprintFP

##  1 Introduction

Symmetries play a central role in our understanding of dynamical processes in
both classical goldstein1980classical ; arnold2013mathematical  and quantum
Coleman:1985rnk  physics. Emmy Noether achieved groundbreaking insight, when
she proved that the presence of a global continuous symmetry in the action
${\cal S}$ of a system implies the existence of a conserved current, whenever
the equations of motions are fulfilled noether1971invariant . Via such a
Noether current, one can define a quantity, which remains unchanged during the
evolution of the system and which is referred to as Noether charge. Noether’s
theorem thus offers a fundamental understanding of central tenets of classical
physics, such as energy and momentum conservation, which it relates to the
invariance of physics under translations in time and space respectively.

In quantum theory, the presence of symmetries limits the type of quantum
fluctuations which may occur Coleman:1985rnk , with measurable consequences
for the spectrum of elementary particles and their bound states. The four
Noether currents associated with space and time translations are
conventionally summarized in a quantity called the energy-momentum tensor
$T^{\mu\nu}(x)$, where $\mu$ and $\nu$ refer to spatial and temporal
components. It offers access to vital properties of a system, one pertinent
example being the energy density profile landau2000classical  of a static
charge distribution via the $\varepsilon(x)=T^{00}(x)$ component or the
corresponding electric field-line configuration via the spatial components
$T_{ij}(x)=F_{i\mu}F^{\mu}_{j}-\frac{1}{4}\delta_{ij}F_{\mu\nu}^{2}$ of the
electromagnetic field $F_{\mu\nu}$, referred to as the Maxwell stress tensor
(see e.g. Yanagihara:2018qqg ).

The simulation of dynamical phenomena in classical and quantum systems is
often performed after discretizing space and time on a finite mesh (for a
discussion of discretization in functional spaces see e.g
cockburn2012discontinuous ). Finite difference schemes, formulated in their
modern summation-by-parts (SBP) form (for reviews see e.g. svard2014review ;
fernandez2014review ; lundquist2014sbp ) offer both conceptual and practical
benefits. The SBP approach in both space and time lundquist2014sbp ;
nordstrom2013summation ; nordstrom2016summation  offers proofs of stability
based on the so-called energy method, which can be extended to high-order
schemes in a straight forward fashion. Not only do SBP operators mimic
integration by parts (IBP) exactly in the discretized setting, but in addition
they constitute a cost effective approximation to differential operators on
many mesh types.

The discretization of space and time in its conventional form, i.e.
considering $x$ and $t$ as independent variables, necessarily affects the
symmetry properties of the system at hand (see e.g. the discussion in
JOHNSON1982147 ). Where the continuum theory e.g. admits translations of any
magnitude, i.e. in particular also infinitesimal ones, the discretized theory
on a space-time mesh with grid spacing $\Delta_{\mu}$ only allows one to shift
space and time by that finite amount. In general this entails that a central
condition of Noether’s theorem, the presence of a continuous symmetry, does
not hold and the corresponding continuum Noether charge fails to remain
constant over time. This is particularly concerning with regards to time
translation symmetry and energy conservation, which are closely related to the
stability of the simulation.

Artificial loss of energy is often considered benign, as it is simply a matter
of loosing accuracy. An artificial increase of energy will, as energy is not
bounded from above, eventually lead to a divergence of the simulated dynamics,
characteristic of an unstable scheme. On the other hand, if energy is
conserved, it puts stringent bounds on the growth of the solution. In the
context of symplectic schemes, which conserve energy on average, one can
relate energy conservation directly to the stability of the numerical scheme
(see e.g. e6656c89-b0e4-3df4-8651-29dbb3d55273  and also
nordstrom2023nonlinear ).

One strategy to retain energy conservation for systems with second order
governing equations is to go over to a Hamiltonian approach, where only space
is discretized, while time remains continuous. One converts the equation of
motion of the Lagrange formalism, which is second order in the time derivative
into a set of two equations of motion of first order, after replacing
velocities with the so-called canonical momentum. After this step, a discrete
phase-space volume preserving time stepping may be implemented (c.f. Verlet-
Størmer verlet_computer_1967 ). This approach crucially hinges on the
availability of a Hamiltonian picture, i.e. whether the canonical momenta can
be defined, which may face difficulties in systems with inherent constraints
or requires the choice of a particular gauge, as in Maxwell’s electrodynamics
dirac_generalized_1950 . Another strategy is to determine whether Noether’s
theorem may be salvaged in the presence of a finite grid spacing
anerot_noethers-type_2020 . One may e.g. consider modifications to the
continuum energy expression, which remain conserved, given a particular choice
of difference approximation. However, as the necessary schemes are not of SBP
type, they do not mimic other relevant properties of the continuum theory.

In this study we develop a generic approach to discretize second order IVPs on
the level of the system Lagrangian, while retaining the manifest translation
invariance of the continuum theory. In order to do so we will take inspiration
from the general theory of relativity (for a textbook see e.g.
stephani2004relativity ), where space and time are treated on the same
footing. In this formalism the presence of translation symmetry is evident
from the form of the Langrangian itself. We build upon our prior work on
formulating IVPs directly via the action of the system, which allows us to
avoid the need to derive their equation of motion. The action of the system is
discretized using SBP finite difference operators with a physical null-space,
developed in our previous paper Rothkopf:2022zfb . These operators are crucial
in mimicking the continuum derivation of Noether’s theorem (and if one wishes
to do so, the equations of motion).

The central outcome of this proof-of-principle study is a prescription of how
to discretize second order IVPs directly on the level of the Lagrangian, while
retaining the continuum time translation symmetry and thus exact conservation
of the corresponding Noether charge. No reference to a Hamiltonian is
required. We observe that a non-equidistant discretization emerges in the time
coordinate, which represents a form of automatic adaptive mesh refinement
(AMR) berger1984adaptive ; lohner1987adaptive ; berger1989local , guided by
the inherent symmetries of the system. Our results open up a novel route for
obtaining optimal AMR procedures, where clustering and coarsening emerge as
part of the solution process, thus avoiding the conventional use of sensors
(see e.g. persson2006sub ), adjoint techniques (see e.g. nemec2008adjoint ;
offermans2023error ) or error estimates (see e.g. mavriplis1994adaptive ;
henderson1999adaptive ; kompenhans2016adaptation ).

In section 2 we discuss the continuum formulation of our geometrized
variational approach with time considered as dependent variable. In section 3
the discretized formalism is introduced and we present its efficacy in section
4 using different example systems. We close with a summary and outlook in
section 5.

##  2 Continuum formalism with manifest translation symmetry

The common starting point for the formulation of the variational principle in
classical point mechanics is to consider the dynamics of a system as boundary
value problem (BVP). The system, which takes on position $x_{i}$ at $t_{i}$
evolves to position $x_{f}$ at $t_{f}$ and we wish to determine the trajectory
it follows. Obviously this formulation is not causal, as we already need to
know the end-point of the dynamics to determine the trajectory. As discussed
in galley_classical_2013  and in our previous study Rothkopf:2022zfb  it is
possible to formulate the variational problem as a genuine initial value
problem through a doubling of the degrees of freedom of the system.

In order to focus on the qualitatively novel ingredients of our variational
approach, we first introduce it in the standard context of point mechanics as
a BVP. The implementation for a genuine IVP is given in the subsequent
subsection.

###  2.1 Boundary value problem formulation

Symmetry is a central mathematical pillar of the theory of relativity. In the
special theory of relativity one formulates the laws of physics in a way that
remains invariant under so-called Lorentz transformations of the coordinates,
while in general relativity one constructs a description, which is invariant
under an even larger class of transformations. Such a theory, invariant under
arbitrary differentiable coordinate transformations, is called
reparametrization invariant.

Reparameterization invariance is achieved by considering both space and time
as dynamical degrees of freedom. In this study we are not interested in
determining the dynamical evolution of space-time itself but will simply
borrow this reparametrization invariant formalism of general relativity for
our purposes of obtaining a symmetry preserving discretization. As our prime
example, we set out to describe the dynamics of a point mass in the presence
of a potential. The first step is to convert this physics question into a
purely geometric problem.

In general relativity, the trajectory of a particle, traveling freely in (a
not necessarily flat) space-time described by the metric tensor $g_{\mu\nu}$,
is given by a path that generalizes the notion of the shortest path on the
corresponding space-time manifold. This path is called a geodesic. While the
particle may move in a $(1+d)$ dimensional space-time with $d$ space and one
time direction, its path traces out a one-dimensional submanifold, which we
can parameterize with a single, so called world-line parameter, denoted in the
following by $\gamma$. We will restrict ourselves here to two dimensions, i.e.
$d=1$, a system with one spatial and one temporal direction expressed in
coordinates as ${\bf x}(\gamma)=(t(\gamma),x(\gamma))$.

A geodesic may be obtained from a variational principle jost1998calculus ,
which asks for the critical point of the following action functional that
measures the length of the path between two space-time points ${\bf
x}(\gamma_{i})$ and ${\bf x}(\gamma_{f})$

$\displaystyle{\cal
S}=\int_{\gamma_{i}}^{\gamma_{f}}\,d\gamma\,(-mc)\sqrt{g_{\mu\nu}\frac{dx^{\mu}}{d\gamma}\frac{dx^{\nu}}{d\gamma}},\qquad{\bf
x}(\gamma_{i})={\bf x}_{i},\;{\bf x}(\gamma_{f})={\bf x}_{f}.$ (1)

Here Einstein’s summation convention has been adopted and we have included the
dimensionful prefactor $mc$, which, as we will show explicitly below, allows
us to recover the usual action in the non-relativistic limit from eq. 1.

We refer to time $t(\gamma)$ as the zeroth component $x^{0}$ of the vector
${\bf x}$ and to the spatial coordinate $x(\gamma)$ as the first component
$x^{1}$. Note that this functional is reparametrization invariant under any
differentiable redefinition of the parameter $\gamma$. I.e. when converting
from $\gamma\to\gamma^{\prime}$ the conversion of differentials under the
square root produces terms $d\gamma^{\prime}/d\gamma$ that cancel with the
conversion factor of the measure.

The geodesics of flat space-time, described by the diagonal metric tensor
$g={\rm diag}[c^{2},-1]$, which arise from the critical point of the action
functional

$\displaystyle{\cal S}_{\rm
flat}=\int_{\gamma_{i}}^{\gamma_{f}}d\gamma(-mc)\sqrt{c^{2}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}-\Big{(}\frac{dx}{d\gamma}\Big{)}^{2}},$
(2)

are straight lines, which are traversed with constant speed (see chapter 3.4
of carroll2019spacetime ), in agreement with Newtonian mechanics.

It is important to note that while our intuition of the concept of shortest
path relies on geometries with positive definite metrics (Riemannian
geometry), physical spacetime, as confirmed by experiment, has a metric with
both positive and negative eigenvalues (pseudo-Riemannian geometry). In such a
geometry the shortest path between two points can denote a saddle point of the
action functional instead of a genuine minimum, as the temporal and spatial
components enter relation (1) with opposite sign.

To describe the presence of an external force acting on a point particle in
flat spacetime, one conventionally amends the action ${\cal S}_{\rm flat}$
simply by adding the potential term $V(x)$ responsible for generating that
force (see chapter 7.9 in goldstein1980classical ).

Let us now discuss how we can exploit the formalism of general relativity to
re-express the evolution of a particle in flat spacetime in the presence of an
external force, instead as an evolution of a free particle in a non-flat
spacetime. In the presence of an external force, encoded in a potential term
$V(x)$, the particle trajectory in flat space-time will deviate from the
straight line. A standard procedure in the study of weak-field gravity is to
reinterpret the change in the particle trajectory due to a potential, instead,
as the effect of a non-flat space-time without a potential present (see e.g.
chapter 8 of carlip2019general ). This reinterpretation is possible, as long
as the values of the potential are smaller than the rest energy $(mc^{2})$ of
the point mass, a condition which is very well fulfilled for the non-
relativistic systems we are interested in solving.

As we will see in the following, one can introduce the effects of a potential
$V(x)$ on a point particle with mass $m$ in the weak-field limit of general
relativity by modifying the temporal component $g_{00}$ of the diagonal metric
tensor

$\displaystyle g_{00}=c^{2}+2V(x)/m,$ (3)

while keeping $g_{11}=-1$. I.e. one endows the metric with a non-trivial
dependence on the spatial coordinate, trading the absence of an explicit
external force for a non-flat spacetime.

Let us now show that such a modification of the metric indeed recovers the
non-relativistic action of a particle in the presence of the potential $V(x)$.
To this end we insert the modified metric eq. 3 into the geodesic action eq.
1:

$\displaystyle{\cal S}=$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}\,d\gamma\,(-mc)\sqrt{g_{00}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}-\Big{(}\frac{dx}{d\gamma}\Big{)}^{2}}$
(4) $\displaystyle\overset{g_{00}>0}{=}$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}\,d\gamma\,(-mc)\sqrt{g_{00}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}}\sqrt{1-\frac{1}{g_{00}}\underbracket{\Big{(}\frac{dx}{d\gamma}\Big{)}^{2}\Big{(}\frac{dt}{d\gamma}\Big{)}^{-2}}_{(dx/dt)^{2}}}$
(5) $\displaystyle\overset{\frac{dx}{dt}^{2}\ll g_{00}\sim c^{2}}{=}$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}\,d\gamma\,\Big{|}\frac{dt}{d\gamma}\Big{|}\,(-mc)\sqrt{g_{00}}\Big{(}1-\frac{1}{2}\frac{1}{g_{00}}\Big{(}\frac{dx}{d\gamma}\Big{)}^{2}\Big{(}\frac{dt}{d\gamma}\Big{)}^{-2}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+{\cal
O}(\frac{1}{g_{00}^{2}}\Big{(}\frac{dx}{dt}\Big{)}^{4})\Big{)}$ (6)
$\displaystyle\overset{V/m\ll c^{2}}{=}$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}\,d\gamma\,\frac{dt}{d\gamma}\Big{(}-mc^{2}+\frac{1}{2}m\Big{(}\frac{dx}{d\gamma}\Big{)}^{2}\Big{(}\frac{dt}{d\gamma}\Big{)}^{-2}-V(x)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+{\cal
O}(\Big{(}\frac{V}{mc^{2}}\Big{)}^{2})+{\cal
O}(\frac{1}{c^{2}}\Big{(}\frac{dx}{dt}\Big{)}^{4})\Big{)}$ (7)
$\displaystyle=$
$\displaystyle\int_{t_{i}}^{t_{f}}\,dt\Big{(}-mc^{2}+\frac{1}{2}m\Big{(}\frac{dx}{dt}\Big{)}^{2}-V(x)\Big{)}.$
(8)

In the third line we have expanded the rightmost square root in eq. 5,
assuming that the square of the physical velocity $(dx/dt)^{2}$ is much
smaller than $g_{00}$, which is to say that the particle velocity $dx/dt$
itself is much smaller than the speed of light $c$. To go from the third to
the fourth line, we have in addition assumed that the potential is much
smaller than the rest energy of the point particle, which allows us to expand
the term $\sqrt{g_{00}}=\sqrt{c^{2}+2V(x)/m}$ in terms of $V(x)/mc^{2}$. We
will look for solutions where time flows forward and thus have dropped the
absolute value around $dt/d\gamma$ at the beginning of the second to last
line. Note that eq. 8 is nothing but the standard non-relativistic action
goldstein1980classical  for a point particle in the presence of an arbitrary
potential term with the rest energy $mc^{2}$ included.

We have thus successfully related the (artificially constructed) fully
geometric description of the particle in a non-flat spacetime in eq. 5 with
the standard description of a particle propagating in flat spacetime in the
presence of an external potential in eq. 8 in the non-relativistic limit.

We see in eq. 8 that time emerges naturally as the independent variable in
which the action integral is formulated. Of course, choosing time as
independent variable hides the inherent reparametrization invariance, which
persists even in the non-relativistic limit in eq. 7. Interestingly it turns
out that eq. 7 is a generalization of the ad-hoc construction of a
reparametrization invariant non-relativistic action, discussed in standard
textbooks on the calculus of variations (see e.g. jost1998calculus ). Equation
7 includes the rest mass term $-mc^{2}(dt/d\gamma)$, which is missing in the
standard derivation and which in the absence of a potential contributes a
dependence on $(dt/d\gamma)$ that plays a role in obtaining a well-defined
critical point for the time degree of freedom.

The reward for our efforts lies in the fact that eq. 4 is manifestly invariant
under the space-time symmetries of our $(1+1)$ dimensional system. If $V(x)=0$
only the derivatives $dt/d\gamma$ and $dx/d\gamma$ but not $t$ and $x$ itself
appear in the action functional eq. 1. In turn adding a constant shift to
either $t$ or $x$ as in ${\bf x}\to{\bf x}+{\bf s}$ leaves the action
invariant. In the presence of a spatially dependent potential $V(x)$,
$g_{00}(x)$ too becomes dependent on space $x$ and only time translation
invariance remains (as the force induced by $V(x)$ changes the momentum of the
point particle).

Proving time translation invariance in the conventional action eq. 8 is much
more involved, as one needs to consider how $x$ as a function of $t$ changes
under such translations and in addition the boundaries of the action integral
themselves are affected by the shift. None of these complications arise in eq.
4111That the derivatives of space and time occur in eq. (4) as squares under
the square root with a relative minus sign (hiding in $g_{11}$) also entails
that the action is manifestly invariant under so called Lorentz boosts. These
transformations mix space and time components and are related to changes
between inertial coordinate systems..

In the calculus of variations it is known that the critical point of the
action ${\cal S}$ can be obtained by solving certain differential equations,
the so called geodesic equations jost1998calculus . It follows from
considering the variation of the action in all of its dependent variables $t$,
$\dot{t}=dt/d\gamma$, $x$ and $\dot{x}=dx/d\gamma$

$\displaystyle\delta{\cal S}[t,\dot{t},x,\dot{x}]=$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\Big{\\{}\frac{\partial{\cal
L}}{\partial t}\delta t+\frac{\partial{\cal
L}}{\partial\dot{t}}\delta\dot{t}+\frac{\partial{\cal L}}{\partial x}\delta
x+\frac{\partial{\cal L}}{\partial\dot{x}}\delta\dot{x}\Big{\\}}$ (9)
$\displaystyle=$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\Big{\\{}\Big{(}\frac{\partial{\cal
L}}{\partial t}-\frac{d}{d\gamma}\frac{\partial{\cal
L}}{\partial\dot{t}}\Big{)}\delta t+\Big{(}\frac{\partial{\cal L}}{\partial
x}-\frac{d}{d\gamma}\frac{\partial{\cal L}}{\partial\dot{x}}\Big{)}\delta
x\Big{\\}}$ (10) $\displaystyle+\left.\Big{[}\frac{\partial{\cal
L}}{\partial\dot{t}}\delta
t\Big{]}\right|_{\gamma_{i}}^{\gamma_{f}}+\left.\Big{[}\frac{\partial{\cal
L}}{\partial\dot{x}}\delta x\Big{]}\right|_{\gamma_{i}}^{\gamma_{f}}.$ (11)

where in the second line we have integrated by parts. As we are considering
the variational problem as boundary value problem with the coordinates $t$ and
$x$ fixed at the start and end points of the trajectory ${\bf
x}(\gamma_{i})={\bf x}_{i},{\bf x}(\gamma_{f})={\bf x}_{f}$, also the
variations $\delta t$ and $\delta x$ on the boundary vanish and so do the two
boundary terms above. Note that we consider $t$ and $x$ as distinct degrees of
freedom, so that the terms in the parentheses, multiplying the arbitrary
variations $\delta x$ and $\delta t$, must vanish each independently at the
stationary point $\delta{\cal S}=0$.

By deriving the Euler-Lagrange equations of the system in the spirit of the
standard BVP treatment of classical mechanics, the above derivation tells us
that we may locate the classical trajectory of a non-relativistic particle
under the influence of a potential, by finding the critical point of the
action eq. 1 with modified $g_{00}$ component of the metric, while keeping the
start and end coordinates ${\bf x}(\gamma_{i})$ and ${\bf x}(\gamma_{f})$
fixed.

Note that there exist infinitely many different parameterizations of the
trajectory described by $\delta{\cal S}=0$, which all differ by the velocity
in $\gamma$, in which this trajectory is traversed. In practice these
different stationary points of ${\cal S}$ lead to difficulties in numerical
optimization and we therefore follow the standard practice (see e.g.
discussion in rizzuti_square_2019  or carroll2019spacetime ) of selecting a
particular parameterization by choosing instead of $\cal S$ the variations of
the functional

$\displaystyle{\cal E}_{\rm BVP}=\int_{\gamma_{i}}^{\gamma_{f}}d\gamma E_{\rm
BVP}[t,\dot{t},x,\dot{x}]=\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\frac{1}{2}\Big{(}g_{00}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}+g_{11}\Big{(}\frac{dx}{d\gamma}\Big{)}^{2}\Big{)}.$
(12)

It differs from $\cal S$ via squaring the integrand and replacing the pre-
factor $-mc$ by $1/2$. These are both irrelevant changes with respect to the
classical equation of motion. Since $\cal E_{\rm BVP}$ and ${\cal S}$ differ
by a monotonous function applied to their integrands, formally the same
critical point ensues. I.e. the variation of $\cal E_{\rm BVP}$ is given by
$\delta{\cal L}=\delta\sqrt{E_{\rm BVP}}=\frac{\delta E_{\rm
BVP}}{2\sqrt{E_{\rm BVP}}}=0$, so that the trajectory that extremizes ${\cal
E}_{\rm BVP}$ agrees with that for ${\cal S}$ at the critical point. Note that
the functional ${\cal E}_{\rm BVP}$ is not reparametrization invariant
anymore. The derivative terms enter quadratically, and produce a conversion
factor $(d\gamma^{\prime}/d\gamma)^{2}$, which cannot be absorbed by the
measure $d\gamma$ alone.

Let us compute the Euler-Lagrange equations (the geodesic equations) for time
$t$ and space $x$ following from the variation of eq. 12

$\displaystyle\delta{\cal E}_{\rm BVP}[t,\dot{t},x,\dot{x}]$ $\displaystyle=$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\Big{\\{}\frac{\partial
E_{\rm BVP}}{\partial t}\delta t+\frac{\partial E_{\rm
BVP}}{\partial\dot{t}}\delta\dot{t}+\frac{\partial E_{\rm BVP}}{\partial
x}\delta x+\frac{\partial E_{\rm BVP}}{\partial\dot{x}}\delta\dot{x}\Big{\\}}$
(13) $\displaystyle=$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\Big{\\{}\Big{(}\frac{\partial
E_{\rm BVP}}{\partial t}-\frac{d}{d\gamma}\frac{\partial E_{\rm
BVP}}{\partial\dot{t}}\Big{)}\delta t+\Big{(}\frac{\partial E_{\rm
BVP}}{\partial x}-\frac{d}{d\gamma}\frac{\partial E_{\rm
BVP}}{\partial\dot{x}}\Big{)}\delta x\Big{\\}}$ (14)
$\displaystyle+\left.\Big{[}\frac{\partial E_{\rm BVP}}{\partial\dot{t}}\delta
t\Big{]}\right|_{\gamma_{i}}^{\gamma_{f}}+\left.\Big{[}\frac{\partial E_{\rm
BVP}}{\partial\dot{x}}\delta x\Big{]}\right|_{\gamma_{i}}^{\gamma_{f}}.$ (15)

As the above boundary terms vanish, we are left with evaluating the individual
expressions appearing in the parentheses of eq. 14. Below we evaluate each of
these terms individually

$\displaystyle\frac{\partial E_{\rm BVP}}{\partial t}=0,\quad\frac{\partial
E_{\rm BVP}}{\partial\dot{t}}=g_{00}(x)\frac{dt}{d\gamma},$ (16)
$\displaystyle\frac{\partial E_{\rm BVP}}{\partial
x}=\frac{1}{2}\frac{\partial g_{00}(x)}{\partial
x}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2},\quad\frac{\partial E_{\rm
BVP}}{\partial\dot{x}}=g_{11}\frac{dx}{d\gamma}=-\frac{dx}{d\gamma},$ (17)

making explicit the ingredients to the geodesic equations for the temporal and
spatial degrees of freedom

$\displaystyle\frac{d}{d\gamma}\Big{(}g_{00}\frac{dt}{d\gamma}\Big{)}=0,$ (18)
$\displaystyle\frac{d}{d\gamma}\Big{(}\frac{dx}{d\gamma}\Big{)}+\frac{1}{2}\frac{\partial
g_{00}}{\partial x}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}=0.$ (19)

The attentive reader will have recognized that eq. 18 constitutes a
conservation equation for the expression inside the parenthesis. In the next
chapter we will show that this quantity indeed is the conserved charge
associated with the time translation symmetry of our system. In general the
geodesic equations do not single out the conserved quantities in such a simple
fashion. There however exists an systematic procedure to identify the space-
time symmetries of the system in the form of different so-called Killing
vectors, each of which leads to one conserved quantity (see section 2.2).

Note that the geodesic equations eqs. 18 and 19 are often written in a more
concise fashion in the general relativity literature (see e.g.
stephani2004relativity ). They are expressed for a general metric using the
so-called Christoffel symbols
$\Gamma^{\alpha}_{\mu\nu}=\frac{1}{2}g^{\alpha\beta}\big{(}\partial
g_{\beta\mu}/\partial x_{\nu}+\partial g_{\beta\nu}/\partial x_{\mu}-\partial
g_{\mu\nu}/\partial x_{\beta})$, where $g^{\alpha\beta}$ refers to the
components of the inverse of the metric $g_{\alpha\beta}$. One obtains in
short hand notation with Einstein summation implied

$\displaystyle\frac{d^{2}x^{\alpha}}{d\gamma^{2}}+\Gamma^{\alpha}_{\mu\nu}\frac{dx^{\mu}}{d\gamma}\frac{dx^{\nu}}{d\gamma}=0.$
(20)

It is important to note that the derivation of the above expression involves
application of the product rule, which in the discrete setting is not valid.
Therefore even though in the continuum eqs. 18 and 19 and eq. 20 are
equivalent, we will work solely with the former, as only integration by parts
(which is exactly mimicked by summation by parts) has been used in their
derivation.

###  2.2 Conserved quantities, Noether’s theorem and stability

Conservation of momentum and energy in general relativity is conceptually more
involved compared to flat space-time, since the comparison of two quantities
at different space-time points becomes a non-trivial operation due to the
effects of a non-flat metric. However there may exist a vector field
$K^{\mu}(x)$ along which transported quantities remain constant. These vector
fields are known as Killing222For completeness we note that a Killing vector
field $K_{\mu}$ is defined as solution to the Killing equation
$\Big{(}\frac{\partial K_{\mu}}{\partial
x^{\nu}}-\Gamma^{\alpha}_{\mu\nu}K_{\alpha}\Big{)}+\Big{(}\frac{\partial
K_{\nu}}{\partial x^{\mu}}-\Gamma^{\alpha}_{\nu\mu}K_{\alpha}\Big{)}=0.$
vector fields $K^{\mu}(x)$. The Killing vector fields are generators of
infinitesimal isometries of the space-time manifold. Moving all points of the
manifold in the direction of the Killing field leaves the manifold unchanged.

As discussed in standard literature on general relativity (see e.g. chapter
3.8 of carroll2019spacetime ), each Killing vector field $K^{\mu}$ can be used
to define a conserved quantity $Q_{K}$ via the expression

$\displaystyle Q_{K}=g_{\alpha\beta}K^{\alpha}\dot{x}^{\beta}.$ (21)

Computing the change of $Q_{K}$ along a geodesic, parameterized by $\gamma$,
one finds from combining eqs. 20 and 21 and the equation that defines the
Killing vector that $dQ_{K}/d\gamma=0$, i.e. it vanishes. We will give an
explicit example of such a conserved quantity below.

More intuitively, one can think of the role of $K^{\mu}$ as pointing out
directions along which the metric $g$ of spacetime in our system remains
constant. In the spirit of Noether’s theorem, assume that the integrand
$E_{\rm BVP}$ of our action functional ${\cal E}_{\rm BVP}$ in eq. 12 remains
unchanged under infinitesimal translations with magnitude $\epsilon$ in the
direction of $K^{\mu}$. The change in coordinates under such a shift is
$\delta x^{\mu}=\epsilon K^{\mu}$. Noether’s theorem tells us that the
conserved quantity corresponding to $\delta x^{\mu}$ is given by $J=\delta
x^{\mu}\frac{\partial E}{\partial\dot{x}^{\mu}}$, which, when written
explicitly as $\epsilon K^{\alpha}g_{\alpha\beta}\frac{dx^{\beta}}{d\gamma}$,
turns out to just be $\epsilon Q_{K}$.

In case of our geometrized problem of determining the dynamics of a point
particle under the influence of a potential $V(x)$, the metric remains
independent of time $t$. Thus the vector $K_{t}=(1,0)$ constitutes a Killing
vector associated with time translation symmetry. The conservation of the
associated conserved quantity
$Q_{t}=K_{t}^{\mu}g_{\mu\nu}\dot{x}^{\nu}=g_{00}\dot{t}$ follows straight
forwardly from the geodesic equation for $t$

$\displaystyle\frac{d}{d\gamma}Q_{t}\overset{K_{t}=(1,0)}{\underset{\lx@cref{creftype~refnum}{eq:killingcons}}{=}}\frac{d}{d\gamma}\Big{(}g_{00}\dot{t}\Big{)}\overset{\lx@cref{creftype~refnum}{eq:geodt}}{=}0,$
(22)

i.e. the quantity $Q_{t}$ remains constant along the geodesic. Note that this
quantity is different from the usual energy considered in the non-relativistic
formalism.

Turning to the question of stability, let us show next that as a consequence
of the presence of a conserved quantity together with the form of the geodesic
equations and the reasonable assumption that the potential of the system is
bounded from below, it is possible to provide an upper bound on the
derivatives of the trajectories obtained as critical point of the functional
eq. 12.

In an analogy to the construction of a Hamiltonian from a Lagrangian, we
define the following

$\displaystyle{\cal H}_{\rm BVP}$
$\displaystyle=\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\underbracket{\frac{1}{2}\Big{(}g_{00}(x)\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}-g_{11}\Big{(}\frac{dx}{d\gamma}\Big{)}^{2}\Big{)}}_{H_{\rm
BVP}},$ (23)
$\displaystyle=\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\frac{1}{2}\Big{(}\big{(}c^{2}+2\frac{V(x)}{m}\big{)}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}+\Big{(}\frac{dx}{d\gamma}\Big{)}^{2}\Big{)}.$
(24)

Due to the flipped sign in front of $g_{11}$, compared to the action eq. 12,
this quantity is actually positive definite, as long as $V(x)$ is bounded from
below333Since physical forces arise from the derivative of the potential, we
may always add a constant to a bounded potential that will make $g_{00}$
positive.. ${\cal H}_{\rm BVP}$ thus provides a norm on the function space in
which $t(\gamma)$ and $x(\gamma)$ reside. Now let us inspect the evolution of
the integrand $H_{\rm BVP}$

$\displaystyle\frac{dH_{\rm BVP}}{d\gamma}$
$\displaystyle=\frac{1}{2}\frac{dg_{00}}{d\gamma}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}+g_{00}\frac{dt}{d\gamma}\frac{d^{2}t}{d\gamma^{2}}+\frac{dx}{d\gamma}\frac{d^{2}x}{d\gamma^{2}},$
$\displaystyle=\frac{dx}{d\gamma}\Big{[}\frac{1}{2}\frac{\partial
g_{00}}{\partial
x}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}+\frac{d^{2}x}{d\gamma^{2}}\Big{]}+g_{00}\frac{dt}{d\gamma}\frac{d^{2}t}{d\gamma^{2}},$
$\displaystyle\overset{\lx@cref{creftype~refnum}{eq:geodx}}{=}\underbracket{g_{00}\frac{dt}{d\gamma}}_{Q_{t}\,\rm
const.}\frac{d^{2}t}{d\gamma^{2}}.$ (25)

To arrive at the final expression in eq. 25, we use the fact that one can
rewrite $dg_{00}/d\gamma=(\partial g_{00}(x)/\partial x)\dot{x}$ and combine
the first and third term to apply eq. 19. This simplification tells us that
the change in $H_{\rm BVP}$ is given solely by the second derivative of time
with respect to the world-line parameter. Now we can integrate up twice ${\cal
H}_{\rm
BVP}=\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\int_{\gamma_{i}}^{\gamma}d\gamma^{\prime}(dH_{\rm
BVP}(\gamma^{\prime})/d\gamma^{\prime})$ to get

$\displaystyle{\cal H}_{\rm BVP}$
$\displaystyle=mg_{00}(x_{i})\dot{t}(\gamma_{i})\;\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\big{(}\dot{t}(\gamma)-\dot{t}(\gamma_{i})\big{)},$
$\displaystyle=g_{00}(x_{i})\dot{t}(\gamma_{i})\;\Big{(}-\dot{t}(\gamma_{i})(\gamma_{f}-\gamma_{i})+(t(\gamma_{f})-t(\gamma_{i}))\Big{)},$
$\displaystyle\leq
g_{00}(x_{i})\dot{t}(\gamma_{i})\big{(}t(\gamma_{f})-t(\gamma_{i})\big{)}.$
(26)

For the last inequality we use the fact that the world-line is parameterized
by an increasing $\gamma$ and correspondingly time moves forward along the
world-line.

In the BVP setting, where both $t(\gamma_{i})$ and $t(\gamma_{f})$ are given
apriori, eq. 26 constitutes a proof that the norm ${\cal H}_{\rm BVP}$ defined
on the derivatives of the solution $t$ and $x$ grows at most linearly with
time, precluding the occurrence of exponentially increasing behavior that
would signal an instability, in turn establishing stability of the geometric
approach.

###  2.3 Initial value formulation

So far we have shown how the geodesic equations eqs. 18 and 19 can be obtained
from a variational principle formulated as a boundary value problem in time.
However for a causal description as an initial value problem, we must be able
to determine the dynamics of the particle without knowledge of the final point
of the trajectory. If one wishes to prescribe only initial values, i.e.
positions and derivatives at $\gamma_{i}$, then the variations $\delta
x^{\mu}$ in eq. 9 do not vanish at the end of the particle world line, i.e. at
$\gamma_{f}$. In turn the equivalence between the critical point of ${\cal S}$
and the Euler-Lagrange equations in eq. 10 does not hold. As discussed by
galley_classical_2013  and put into practice in our previous publication
Rothkopf:2022zfb  one can overcome this issue by constructing an action with
doubled degrees of freedom, living on a closed contour with a forward and
backward branch in $\gamma$.

Since both time and space constitute dependent degrees of freedom in our
approach, we need to introduce both forward and backward variants of each of
them $x_{1}(\gamma),x_{2}(\gamma)$ and $t_{1}(\gamma),t_{2}(\gamma)$. The
degrees of freedom on the forward contour enter the action functional with the
usual Lagrangian, while those on the backward contour are assigned the
negative Lagrangian. Choosing to build the doubled formalism based on the
action ${\cal E}_{\rm BVP}$ we obtain

$\displaystyle{\cal E}_{\rm IVP}$
$\displaystyle=\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\,E_{\rm
IVP}[t_{1},\dot{t}_{1},x_{1},\dot{x}_{1},t_{2},\dot{t}_{2},x_{2},\dot{x}_{2}],$
(27) $\displaystyle=\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\Big{\\{}E_{\rm
BVP}[t_{1},\dot{t}_{1},x_{1},\dot{x}_{1}]-E_{\rm
BVP}[t_{2},\dot{t}_{2},x_{2},\dot{x}_{2}]\Big{\\}}.$ (28)

As discussed in detail in Rothkopf:2022zfb , the inner workings of the doubled
formalism become more transparent, once we go over to expressing the action
$\cal E_{\rm IVP}$ in terms of the central and difference coordinates
$x_{+}=\frac{1}{2}(x_{1}+x_{2})$ and $x_{-}=x_{1}-x_{2}$ and
$t_{+}=\frac{1}{2}(t_{1}+t_{2})$ and $t_{-}=t_{1}-t_{2}$ respectively. The
variation now proceeds in the independent degrees of freedom $x_{\pm}$ and
$t_{\pm}$ and yields

$\displaystyle\delta{\cal E}_{\rm
IVP}[t_{\pm},\dot{t}_{\pm},x_{\pm},\dot{x}_{\pm}]=$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\Big{\\{}\frac{\partial
E_{\rm IVP}}{\partial t_{+}}\delta t_{+}+\frac{\partial E_{\rm
IVP}}{\partial\dot{t}_{+}}\delta\dot{t}_{+}+\frac{\partial E_{\rm
IVP}}{\partial t_{-}}\delta t_{-}+\frac{\partial E_{\rm
IVP}}{\partial\dot{t}_{-}}\delta\dot{t}_{-}$ (29)
$\displaystyle\qquad+\frac{\partial E_{\rm IVP}}{\partial x_{+}}\delta
x_{+}+\frac{\partial E_{\rm
IVP}}{\partial\dot{x}_{+}}\delta\dot{x}_{+}+\frac{\partial E_{\rm
IVP}}{\partial x_{-}}\delta x_{-}+\frac{\partial E_{\rm
IVP}}{\partial\dot{x}_{-}}\delta\dot{x}_{-}\Big{\\}}$ $\displaystyle=$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\Big{\\{}\Big{(}\frac{\partial
E_{\rm IVP}}{\partial t_{+}}-\frac{d}{d\gamma}\frac{\partial E_{\rm
IVP}}{\partial\dot{t}_{+}}\Big{)}\delta t_{+}+\Big{(}\frac{\partial E_{\rm
IVP}}{\partial t_{-}}-\frac{d}{d\gamma}\frac{\partial E_{\rm
IVP}}{\partial\dot{t}_{-}}\Big{)}\delta t_{-}$
$\displaystyle+\Big{(}\frac{\partial E_{\rm IVP}}{\partial
x_{+}}-\frac{d}{d\gamma}\frac{\partial E_{\rm
IVP}}{\partial\dot{x}_{+}}\Big{)}\delta x_{+}+\Big{(}\frac{\partial E_{\rm
IVP}}{\partial x_{-}}-\frac{d}{d\gamma}\frac{\partial E_{\rm
IVP}}{\partial\dot{x}_{-}}\Big{)}\delta x_{-}\Big{\\}}$ (30)
$\displaystyle+\left.\Big{[}\frac{\partial E_{\rm
IVP}}{\partial\dot{t}_{+}}\delta
t_{+}\Big{]}\right|_{\gamma_{i}}^{\gamma_{f}}+\left.\Big{[}\frac{\partial
E_{\rm IVP}}{\partial\dot{t}_{-}}\delta
t_{-}\Big{]}\right|_{\gamma_{i}}^{\gamma_{f}}+\left.\Big{[}\frac{\partial
E_{\rm IVP}}{\partial\dot{x}_{+}}\delta
x_{+}\Big{]}\right|_{\gamma_{i}}^{\gamma_{f}}+\left.\Big{[}\frac{\partial
E_{\rm IVP}}{\partial\dot{x}_{-}}\delta
x_{-}\Big{]}\right|_{\gamma_{i}}^{\gamma_{f}}.$

To arrive at eq. 30 we have carried out four integrations by parts. As the
next step, we consider under which conditions the boundary terms in the above
expression vanish. Since we prescribe fixed initial values for both time and
space, the variations $\delta t_{\pm}(\gamma_{i})=0$ and $\delta
x_{\pm}(\gamma_{i})=0$ vanish. What about the variations at the end of the
forward and backward world-line? As long as we require that

$\displaystyle x_{2}(\gamma_{f})=x_{1}(\gamma_{f}),\quad
t_{2}(\gamma_{f})=t_{1}(\gamma_{f}),$ (31)

it follows that $\delta x_{-}(\gamma_{f})$ and $\delta t_{-}(\gamma_{f})$
vanish and with it the corresponding boundary terms. The only remaining terms
are those at $\gamma_{f}$ which feature $\delta x_{+}$ and $\delta t_{+}$. As
these variations do not vanish, we instead inspect the terms multiplying them,
i.e. $\partial E_{\rm IVP}/\partial\dot{t}_{+}$ and $\partial E_{\rm
IVP}/\partial\dot{x}_{+}$. Using the definition $x_{1}=x_{+}+\frac{1}{2}x_{-}$
and $x_{2}=x_{+}-\frac{1}{2}x_{-}$ and correspondingly for $t_{1,2}$, we find
from the defining equation for $E_{\rm IVP}$ eq. 12

$\displaystyle\frac{dE_{\rm IVP}}{d\dot{x}_{+}}$ $\displaystyle=\frac{\partial
E_{\rm
IVP}[t_{1,2},\dot{t}_{1,2},x_{1,2},\dot{x}_{1,2}]}{\partial\dot{x}_{1}}\frac{d\dot{x}_{1}}{d\dot{x}_{+}}+\frac{\partial
E_{\rm
IVP}[t_{1,2},\dot{t}_{1,2},x_{1,2},\dot{x}_{1,2}]}{\partial\dot{x}_{2}}\frac{d\dot{x}_{2}}{d\dot{x}_{+}},$
$\displaystyle=\frac{\partial E_{\rm
BVP}[t_{1},\dot{t}_{1},x_{1},\dot{x}_{1}]}{\partial\dot{x}_{1}}\frac{d\dot{x}_{1}}{d\dot{x}_{+}}-\frac{\partial
E_{\rm
BVP}[t_{2},\dot{t}_{2},x_{2},\dot{x}_{2}]}{\partial\dot{x}_{2}}\frac{d\dot{x}_{2}}{d\dot{x}_{+}},$
(32)
$\displaystyle=g_{11}(x_{1})\dot{x}_{1}-g_{11}(x_{2})\dot{x}_{2}=-\dot{x}_{1}+\dot{x}_{2}.$
(33)

Similarly one obtains

$\displaystyle\frac{dE_{\rm
IVP}}{d\dot{t}_{+}}=g_{00}(x_{1})\dot{t}_{1}-g_{00}(x_{2})\dot{t}_{2}.$ (34)

Together with condition eq. 31 that the values of $x_{1,2}$ and $t_{1,2}$ must
agree at $\gamma_{f}$, this result tells us that in order for the two
remaining boundary terms to vanish, we need to also identify the derivatives
of $x_{1,2}$ and $t_{1,2}$ at the point $\gamma_{f}$

$\displaystyle\dot{x}_{2}(\gamma_{f})=\dot{x}_{1}(\gamma_{f}),\quad\dot{t}_{2}(\gamma_{f})=\dot{t}_{1}(\gamma_{f}).$
(35)

Note that we have now managed to remove the boundary terms without the need
for specifying the concrete value of t’s and x’s at the final point
$\gamma_{f}$. This is the central contribution of the forward-backward
construction.

The last remaining step is to undo the proliferation of degrees of freedom
that occurred when introducing the forward-backward construction. It has been
shown Berges:2007ym ; galley_classical_2013  that taking the so-called
physical limit achieves this goal, where the constraints
$x_{1}(\gamma)-x_{2}(\gamma)=x_{-}(\gamma)=0$ and
$t_{1}(\gamma)-t_{2}(\gamma)=t_{-}(\gamma)=0$ are enforced. The remaining
$x_{+}$ and $t_{+}$ are identified with the true classical geodesics.

In terms of the Euler-Lagrange equations in parentheses in eq. 30

$\displaystyle\frac{\partial E_{\rm IVP}}{\partial
x_{\pm}}-\frac{d}{d\gamma}\frac{\partial E_{\rm
IVP}}{\partial\dot{x}_{\pm}}=0,\qquad\frac{\partial E_{\rm IVP}}{\partial
t_{\pm}}-\frac{d}{d\gamma}\frac{\partial E_{\rm
IVP}}{\partial\dot{t}_{\pm}}=0,$ (36)

the physical limit entails that only those equations independent of $x_{-}$
and $t_{-}$ survive. With the construction of the action $E_{\rm IVP}=E_{\rm
BVP}[x_{1},\dot{x}_{1},t_{1},\dot{t}_{1}]-E_{\rm
BVP}[x_{2},\dot{x}_{2},t_{2},\dot{t}_{2}]$ from a difference of the $E_{\rm
BVP}$ functionals, there will appear at least a linear dependence on the minus
degrees of freedom. Hence in the physical limit only those Euler-Lagrange
equations linear in $x_{-}$ and $t_{-}$ will survive, where the minus degrees
of freedom have been removed by taking the derivative with respect to $x_{-}$
or $t_{-}$.

Note that we have decided to not only specify the value and derivative of $x$
at initial $\gamma_{i}$ but also those of $t$. As we wish to determine the
dynamics of a point particle in the presence of a potential with given
$x(t_{i})$ and $dx/dt(t_{i})$, there remains a freedom in choosing
$\dot{x}(\gamma_{i})$ and $\dot{t}(\gamma_{i})$, since only their ratio needs
to be fixed $dx/dt(t=t_{0})=\dot{x}(\gamma_{i})/\dot{t}(\gamma_{i})$. The end
of the time interval traversed by the world line parameter $\gamma$, will
consequently depend on the value prescribed to $\dot{t}(\gamma_{i})$ and
emerges dynamically from the combined evolution of $x$ and $t$.

At this point we have formulated a manifest time translation symmetric
variational principle that encodes the dynamics of a point particle evolving
in the presence of a non-relativistic potential as initial value problem. Our
next goal is to discretize the action functional ${\cal E}_{\rm IVP}$ in
section 3 using SBP finite difference operators. Since all derivations of the
Euler-Lagrange equations, as well as that of the conserved quantity $Q_{t}$
have made ample reference to integration by parts, it is paramount to use such
a discretization technique, which faithfully mimics this continuum property on
a finite mesh.

##  3 Discretized formalism for IVPs

The central novelty we introduce in this section is related to the fact that
the discretization of the action functional takes place in the world-line
parameter $\gamma$ and not in the time variable $t$, as in conventional
discretization prescriptions. I.e. the values of both time $t(\gamma)$ and
position $x(\gamma)$ remain continuous and in turn we achieve preservation of
the continuum space-time symmetries even after discretization.

In the presence of a potential that depends on $x$ but not on $t$, the
invariance under infinitesimal constant shifts in time is hence retained. This
comes about, since the metric remains invariant under changes in $t$, which in
turn leads to a simple form of the corresponding Killing equation, which shows
that $K_{t}=(1,0)$ indeed is a Killing vector. The symmetry of the metric
under time translation is intimately related to energy conservation via
$Q_{t}$ and thus the stability of the simulation. In the absence of a
potential, when the metric does not depend on neither $t$ nor $x$, our
discretized approach, in addition to $K_{t}=(1,0)$, retains the continuum
invariance under shifts in $x$ via the Killing vector $K_{x}=(0,1)$, as well
as the invariance under boosts via the Killing vector $K_{\eta}=(x,t)$.

We will give numerical evidence that we achieve exact conservation of $Q_{t}$
in the interior of the simulated domain, even in the case of highly non-
harmonic motion. In contrast to other formally energy preserving schemes, such
as the leap-frog, our approach, using SBP operators, is consistent with the
continuum formulation, in that it only requires the actual initial conditions
of the system at hand, avoiding the need to stagger the degrees of freedom
(also known as insertion of dummy points).

After introducing the discretization on the level of the underlying action
functional, we will obtain the classical trajectory by numerically finding the
critical point of that functional without the need to derive the corresponding
equations of motion. To make sure that the solution of the discretized
variational principle mimics as accurately as possible the continuum theory,
we deploy summation-by-parts finite difference operators svard2014review ;
fernandez2014review ; lundquist2014sbp .

Note that we are discretizing the world-line parameter $\gamma$ with
equidistant steps, whereas both the values of $t$ and $x$ arise dynamically
from the evolution of the simulation along $\gamma$. I.e. a not necessarily
equidistant discretization of the time coordinate emerges dynamically in our
approach. As we will see in section 4 this dynamical time discretization
realizes a one-dimensional form of automatic adaptive mesh refinement, guided
by the symmetries of the system. I.e. the non-equidistant discretization in
$t$ plays a crucial role in guaranteeing that the Noether charge $Q_{t}$
remains conserved.

Another non-standard feature of our technique is the departure from the
conventional notion of carrying out a simulation on a predefined time
interval. We instead provide the initial time and its velocity with respect to
$\gamma$, so that the end-point of the simulation too emerges dynamically.

In the following we will consider the trajectory of a point particle
propagating under the influence of an arbitrary $x$ but not $t$ dependent
potential $V(x)$. We begin by discretizing the action functional ${\cal
E}_{\rm IVP}$ of eq. 28 along the world-line parameter $\gamma$ between
$\gamma_{i}$ and $\gamma_{f}$ with $N_{\gamma}$ steps, leading to a step-size
of $d\gamma=(\gamma_{f}-\gamma_{i})/(N_{\gamma}-1)$. We will add to $E_{\rm
IVP}$ Lagrange multipliers to explicitly account for both the initial
conditions and the connecting conditions required by doubling of the degrees
of freedom. The forward and backward paths $x_{1,2}$ and times $t_{1,2}$ are
described by ${\bf
x}_{1,2}=(x_{1,2}(0),x_{1,2}(\Delta\gamma),x_{1,2}(2\Delta\gamma),\ldots,x_{1,2}((N_{\gamma}-1)\Delta\Gamma))^{\rm
T}$ and ${\bf
t}_{1,2}=(t_{1,2}(0),t_{1,2}(\Delta\gamma),t_{1,2}(2\Delta\gamma),\ldots,t_{1,2}((N_{\gamma}-1)\Delta\gamma))^{\rm
T}$ respectively.

The integral in ${\cal E}_{\rm IVP}$ is approximated with a quadrature rule,
consistent with our choice of finite difference operator, in the form of a
diagonal positive definite matrix $\mathds{H}$. The inner product on
discretized paths and times thus reads $({\bf x},{\bf x}^{\prime})={\bf
x}^{\rm T}\mathds{H}{\bf x}^{\prime}$.

With integration by parts being a central element in establishing both
equations of motion and the existence of conserved quantities, we must use a
discretization that mimics IBP exactly, which is achieved by deploying
summation-by-parts (SBP) operators $\mathds{D}$ with the defining properties

$\displaystyle\mathds{D}=\mathds{H}^{-1}\mathds{Q},\qquad\mathds{Q}^{\rm
T}+\mathds{Q}=\mathds{E}_{N}-\mathds{E}_{0}={\rm diag}[-1,0,\ldots,0,1].$ (37)

In this study we consider both the lowest order SBP discretization scheme,
referred to as SBP21 and the next higher order scheme SBP42. The former is
second order in the interior and exhibits one order less on the boundary.
Using the trapezoidal rule for integration one has

$\mathds{H}^{[2,1]}=\Delta\gamma\left[\begin{array}[]{ccccc}1/2&&&&\\\
&1&&&\\\ &&\ddots&&\\\ &&&1&\\\
&&&&1/2\end{array}\right],\quad\mathds{D}^{[2,1]}=\frac{1}{2\Delta\gamma}\left[\begin{array}[]{ccccc}-2&2&&&\\\
-1&0&1&&\\\ &&\ddots&&\\\ &&-1&0&1\\\ &&&-2&2\end{array}\right].$ (38)

The SBP42 scheme achieves fourth order accuracy in the interior, which reduces
to second order on the boundary

$\displaystyle\mathds{H}^{[4,2]}=\Delta\gamma\left[\begin{array}[]{ccccccc}\frac{17}{48}&&&&&&\\\
&\frac{59}{48}&&&&&\\\ &&\frac{43}{48}&&&&\\\ &&&\frac{49}{48}&&&\\\
&&&&&1&\\\ &&&&&&\ddots\\\ \end{array}\right],$ (45)

$\displaystyle\mathds{D}^{[4,2]}=\frac{1}{\Delta\gamma}\left[\begin{array}[]{ccccccccc}-\frac{24}{17}&\frac{59}{34}&-\frac{4}{17}&-\frac{3}{34}&&&&&\\\
-\frac{1}{2}&0&\frac{1}{2}&0&&&&&\\\
\frac{4}{43}&-\frac{59}{86}&0&\frac{59}{86}&-\frac{4}{43}&&&&\\\
\frac{3}{98}&0&-\frac{59}{86}&0&\frac{32}{49}&-\frac{4}{49}&&&\\\
&&\frac{1}{12}&-\frac{2}{3}&0&\frac{2}{3}&-\frac{1}{12}&&\\\
&&&&&&&\ddots&\end{array}\right].$ (52)

The SBP operators defined above are not yet ready for duty in our variational
approach, as they allow for non-physical zero modes. As discussed in detail in
Rothkopf:2022zfb , we can construct null-space consistent444Note that in the
context of PDE’s, SBP operators are considered null-space consistent by
construction, as only their right eigenvectors play a role in the equation of
motion. Here due to the presence of $\mathds{D}^{\rm T}$ in the action
functional, also the left eigenvectors contribute, among which a highly
oscillating null-mode (the so-called $\pi$-mode) can be identified (see ref.
Rothkopf:2022zfb ) SBP operators $\bar{\mathds{D}}$ from the conventional
$\mathds{D}$ by deploying affine coordinates and by absorbing penalty terms,
inspired by the simultaneous-approximation terms (SAT) technique
carpenter1994time , used to regularize SBP operators. A brief overview of this
regularization is given in appendix A.

The idea behind the penalty term construction is that we are assigning a
penalty to all functions that do not fulfill the initial conditions in $t$ and
$x$, which includes the non-physical zero mode of $\mathds{D}$. In turn, when
we will be searching for the critical point of the discretized action
functional ${\cal E}_{\rm IVP}$ the minimizer will approach the correct
solution globally and the presence of the penalty term effectively prevents
contamination of the correct solution by the non-constant zero mode.

Explicitly our regularized and null-space consistent SBP21 operators read

$\displaystyle{\bar{\mathds{D}}}^{{\rm
R},[2,1]}_{t}=\left[\begin{array}[]{cccccc}-\frac{1}{\Delta\gamma}+\frac{2}{\Delta\gamma}&\frac{1}{\Delta\gamma}&&&&-\frac{2}{\Delta\gamma}t_{i}\\\
-\frac{1}{2\Delta\gamma}&0&\frac{1}{2\Delta\gamma}&&&0\\\ &&\ddots&&&\vdots\\\
&&-\frac{1}{2\Delta\gamma}&0&\frac{1}{2\Delta\gamma}&0\\\
&&&-\frac{1}{\Delta\gamma}&\frac{1}{\Delta\gamma}&0\\\ 0&&\ldots&&0&1\\\
\end{array}\right]$ (59) $\displaystyle{\bar{\mathds{D}}}^{{\rm
R},[2,1]}_{x}=\left[\begin{array}[]{cccccc}-\frac{1}{\Delta\gamma}+\frac{2}{\Delta\gamma}&\frac{1}{\Delta\gamma}&&&&-\frac{2}{\Delta\gamma}x_{i}\\\
-\frac{1}{2\Delta\gamma}&0&\frac{1}{2\Delta\gamma}&&&0\\\ &&\ddots&&&\vdots\\\
&&-\frac{1}{2\Delta\gamma}&0&\frac{1}{2\Delta\gamma}&0\\\
&&&-\frac{1}{\Delta\gamma}&\frac{1}{\Delta\gamma}&0\\\ 0&&\ldots&&0&1\\\
\end{array}\right].$ (66)

Using the operators defined above, we can now write the discretized action
functional in the following fashion

$\displaystyle\mathds{E}_{\rm IVP}=$
$\displaystyle\frac{1}{2}\left\\{({\bar{\mathds{D}}^{\rm R}}_{t}{\bf
t}_{1})^{\rm T}\mathbb{d}\left[c^{2}+\frac{2{\bf V}({\bf
x}_{1})}{m}\right]{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{t}{\bf
t}_{1})-({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{1})^{\rm
T}{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{1})\right\\}$
$\displaystyle-$ $\displaystyle\frac{1}{2}\left\\{({\bar{\mathds{D}}^{\rm
R}}_{t}{\bf t}_{2})^{\rm T}\mathbb{d}\left[c^{2}+\frac{2{\bf V}({\bf
x}_{2})}{m}\right]{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{t}{\bf
t}_{2})-({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{2})^{\rm
T}{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{2})\right\\}$
$\displaystyle+$ $\displaystyle\lambda_{1}\big{(}{\bf
t}_{1}[1]-t_{i}\big{)}+\lambda_{2}\big{(}(\mathds{D}{\bf
t}_{1})[1]-\dot{t}_{i}\big{)}+\lambda_{3}\big{(}{\bf x}_{1}[1]-x_{i}\big{)}$
$\displaystyle+$ $\displaystyle\lambda_{4}\big{(}(\mathds{D}{\bf
x}_{1})[1]-\dot{x}_{i}\big{)}$ $\displaystyle+$
$\displaystyle\lambda_{5}\big{(}{\bf t}_{1}[N_{\gamma}]-{\bf
t}_{2}[N_{\gamma}]\big{)}+\lambda_{6}\big{(}{\bf x}_{1}[N_{\gamma}]-{\bf
x}_{2}[N_{\gamma}]\big{)}$ $\displaystyle+$
$\displaystyle\lambda_{7}\big{(}(\mathds{D}{\bf
t}_{1})[N_{\gamma}]-(\mathds{D}{\bf
t}_{2})[N_{\gamma}]\big{)}+\lambda_{8}\big{(}(\mathds{D}{\bf
x}_{1})[N_{\gamma}]-(\mathds{D}{\bf x}_{2})[N_{\gamma}]\big{)}.$ (67)

Conventional matrix vector multiplication is implied in the above expression,
whenever a matrix quantity such as $\bar{\mathds{H}}$ or $\bar{\mathds{D}}$
acts on a vector ${\bf x}_{1,2}$ or ${\bf t}_{1,2}$. The matrix denoted by
$\mathbb{d}[f({\bf x})]$ contains on its diagonal the values
$\mathbb{d}_{kk}=f({\bf x}(\gamma_{k}))$ and zero otherwise. We deploy an
appropriately modified matrix ${\bar{\mathds{H}}}$ for the inner product in
the presence of the affine-coordinate regularized SBP operators (see appendix
A).

The initial conditions we supply are the values of the spatial and temporal
coordinate $x_{i},t_{i}$, as well as the initial velocities with respect to
the world line parameter $\gamma$, i.e. $\dot{x}_{i}$ and $\dot{t}_{i}$. Since
our physical problem is formulated as an initial value problem, given $t_{i}$,
$x_{i}$ and the physical velocity $v_{i}=dx/dt$, there exists a freedom to
choose $\dot{t}_{i}$ and $\dot{x}_{i}$, as only their ratio is fixed
$v_{i}=\dot{x}_{i}/\dot{t}_{i}$. We have added eight Lagrange multipliers,
whose role is to explicitly implement the initial conditions $(\lambda_{1-4})$
and the connecting conditions at the end of the forward and backward branches
of our doubled degree of freedom construction $(\lambda_{5-8})$.

Once the action functional has been formulated in its discrete form, changing
from SBP21 to SBP42 only requires replacement of the corresponding difference
operator $\mathds{D}$ and quadrature matrix $\mathds{H}$ but no further
changes to the functional itself.

This concludes the description of our novel variational approach and we
proceed to evaluate its properties and performance based on two concrete
numerical examples.

##  4 Numerical results

In this section we will present explicit results for the numerically obtained
classical trajectory of a point particle in the presence of two different
potentials, $V_{1}(x)=\alpha x$ and $V_{4}(x)=\kappa x^{4}$. These two choices
correspond to a model of a point mass falling in a constant gravitational
field and carrying out highly-nonlinear anharmonic motion. We set the mass of
the particle to unity, as well as adopt without loss of generality the
convention that the speed of light $c=1$, which simply amounts to a particular
choice of units for length and time.

Let us stress again that while standard numerical methods exist to solve the
equations of motion for each of these systems, the novelty of the approach
presented here lies in the fact that we retain the continuum time shift
invariance of the system and thus achieve exact conservation of $Q_{t}$ in the
interior of the simulated time domain. In addition we determine the classical
trajectory directly from the action functional of the geometrized problem,
without the need to derive the equation of motion.

We implement the action functional eq. 67 in the Mathematica language555The
code using both the SBP21 or SBP42 operator is available under open access on
the Zenodo repository zenodoIVP2023 .. As the critical point of the action may
be a saddle point, instead of an actual minimum, we must be careful in
deploying established numerical optimization algorithms in the dynamical
degrees of freedom ${\bf d}=\\{{\bf t}_{1,2},{\bf x}_{1,2},\lambda_{1-8}\\}$.
Instead of minimizing $\mathds{E}_{\rm IVP}$ directly, we will minimize the
Euclidean norm of the gradient $|\nabla_{\bf d}\mathds{E}_{\rm IVP}|^{2}$. Via
this detour, a saddle point is converted into a minimum. In practice we deploy
a chain of minimization algorithms. We start with a preconditioning based on
the LBFGS quasi-Newton algorithm, which features cost efficient iteration
steps, when far away from the true critical point. It is followed by further
iterations based on the full Newton method, which exhibits a faster
convergence rate than the LBFGS algorithm when close to the critical point.
Once the critical point has been approached to at least floating point
precision we switch to the interior point optimization, which showed reliable
performance in identifying the critical point to any desired tolerance. For
our numerical tests in Mathematica, we used WorkingPrecision of $40$ and
PrecisionGoal of $40$.

The figures shown in the following are based on results from the SBP21
operator and include the outcomes from the SBP42 operators when indicated in
the text.

###  4.1 Linear potential case

We discretize the continuous action functional

$\displaystyle{\cal E}^{\rm lin}_{\rm IVP}=$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\frac{1}{2}\left\\{\Big{(}1+2\alpha
x_{1}(\gamma)\Big{)}\Big{(}\frac{dt_{1}}{d\gamma}\Big{)}^{2}-\Big{(}\frac{dx_{1}}{d\gamma}\Big{)}^{2}\right\\}$
$\displaystyle-$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\frac{1}{2}\left\\{\Big{(}1+2\alpha
x_{2}(\gamma)\Big{)}\Big{(}\frac{dt_{2}}{d\gamma}\Big{)}^{2}-\Big{(}\frac{dx_{2}}{d\gamma}\Big{)}^{2}\right\\}$
$\displaystyle+$
$\displaystyle\lambda_{1}\big{(}t(\gamma_{i})-t_{i}\big{)}+\lambda_{2}\big{(}\dot{t}_{1}(\gamma_{i})-\dot{t}_{i}\big{)}+\lambda_{3}\big{(}x_{1}(\gamma_{i})-x_{i}\big{)}+\lambda_{4}\big{(}\dot{x}_{1}(\gamma_{i})-\dot{x}_{i}\big{)}$
$\displaystyle+$
$\displaystyle\lambda_{5}\big{(}t_{1}(\gamma_{f})-t_{2}(\gamma_{2})\big{)}+\lambda_{6}\big{(}\dot{t}_{1}(\gamma_{f})-\dot{t}_{2}(\gamma_{2})\big{)}$
$\displaystyle+$
$\displaystyle\lambda_{7}\big{(}x_{1}(\gamma_{f})-x_{2}(\gamma_{2})\big{)}+\lambda_{8}\big{(}\dot{x}_{1}(\gamma_{f})-\dot{x}_{2}(\gamma_{2})\big{)}$
(68)

along the world-line of the particle motion between $\gamma_{i}=0$ and
$\gamma_{f}=1$ with $N_{\gamma}=32$ points. Without loss of generality, we
arbitrarily set the starting time to $t_{i}=0$ and the starting position to
$x_{i}=1$. To obtain an initial velocity $v_{i}=1/10$ we choose $\dot{t}=1$
and $\dot{x}=v_{i}$. Note that we do not fix the value of $t_{f}$ but only the
initial velocity of time with respect to $\gamma$. The choice of $\dot{t}=1$
will lead to dynamics, such that $t_{f}$ will be of the order of one. (In the
next subsection we will also provide results for different choices of
$\dot{t}_{i}$.) As strength for the linear potential we choose $\alpha=1/4$.
The corresponding discrete action functional reads explicitly

$\displaystyle\mathds{E}^{\rm lin}_{\rm IVP}=$
$\displaystyle\frac{1}{2}\left\\{({\bar{\mathds{D}}^{\rm R}}_{t}{\bf
t}_{1})^{\rm T}\mathbb{d}\left[1+2\alpha{\bf
x}_{1}\right]{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{t}{\bf
t}_{1})-({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{1})^{\rm
T}{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{1})\right\\}$
$\displaystyle-$ $\displaystyle\frac{1}{2}\left\\{({\bar{\mathds{D}}^{\rm
R}}_{t}{\bf t}_{2})^{\rm T}\mathbb{d}\left[1+2\alpha{\bf
x}_{2}\right]{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{t}{\bf
t}_{2})-({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{2})^{\rm
T}{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{2})\right\\}$
$\displaystyle+$ $\displaystyle\lambda_{1}\big{(}{\bf
t}_{1}[1]-t_{i}\big{)}+\lambda_{2}\big{(}(\mathds{D}{\bf
t}_{1})[1]-\dot{t}_{i}\big{)}$ $\displaystyle+$
$\displaystyle\lambda_{3}\big{(}{\bf
x}_{1}[1]-x_{i}\big{)}+\lambda_{4}\big{(}(\mathds{D}{\bf
x}_{1})[1]-\dot{x}_{i}\big{)}$ $\displaystyle+$
$\displaystyle\lambda_{5}\big{(}{\bf t}_{1}[N_{\gamma}]-{\bf
t}_{2}[N_{\gamma}]\big{)}+\lambda_{6}\big{(}{\bf x}_{1}[N_{\gamma}]-{\bf
x}_{2}[N_{\gamma}]\big{)}$ $\displaystyle+$
$\displaystyle\lambda_{7}\big{(}(\mathds{D}{\bf
t}_{1})[N_{\gamma}]-(\mathds{D}{\bf
t}_{2})[N_{\gamma}]\big{)}+\lambda_{8}\big{(}(\mathds{D}{\bf
x}_{1})[N_{\gamma}]-(\mathds{D}{\bf x}_{2})[N_{\gamma}]\big{)}.$ (69)

Let us take a look in fig. 1 at the raw results for the forward and backward
time and spatial coordinates, as obtained from the critical point of
$\mathds{E}^{\rm lin}_{\rm IVP}$ with $V(x)=\alpha x$. In the top panel, we
show ${\bf t}_{1}(\gamma_{i})$ as red circles and ${\bf t}_{2}(\gamma_{i})$ as
blue crosses, while in the bottom panel these symbols denote the spatial
coordinate of the point particle trajectory ${\bf x}_{1}(\gamma_{i})$ and
${\bf x}_{2}(\gamma_{i})$ respectively. As required by the physical limit
(discussed in section 2.3), we find that the values of the doubled degrees of
freedom coincide at the critical point. The solution of the corresponding
continuum geodesic equations, obtained via the LSODA algorithm of
Mathematica’s NDSolve command is shown as gray solid line and excellent
agreement is observed.

Figure 1: The values of the (top) time coordinates ${\bf t}_{1}(\gamma_{i})$
(red circles) and ${\bf t}_{2}(\gamma_{i})$ (blue crosses) and that of the
(bottom) spatial coordinates ${\bf x}_{1}(\gamma_{i})$ (red circles) and ${\bf
x}_{2}(\gamma_{i})$ (blue crosses) along the world-line parameter $\gamma$, as
obtained from the critical point of $\mathds{E}^{\rm lin}_{\rm IVP}$ with
$V(x)=\alpha x$, discretized with $N_{\gamma}=32$ points and the SBP21
operators. The solution of the corresponding geodesic equations via
Mathematica’s NDSolve is shown as light gray solid line.

Note that due to our choice of $\dot{t}_{i}=1$ the maximum time traversed by
the simulation is close to one.

At first sight it appears that an equidistant discretization of time in
$\gamma$ emerges, but an inspection of the velocity of time with respect to
$\gamma$ in fig. 2 reveals that the time spacing dynamically adapts to the
behavior observed in the spatial coordinate $x$. Close to the maximum of
$x(\gamma)$ at around $\gamma=0.4$ the temporal spacing e.g. has a minimum.
This dynamically emerging time discretization constitutes an automatically
generated non-trivial mesh for the time coordinate and arises naturally in our
formalism. In fact an automatic AMR procedure results.

Figure 2: The change of time with respect to the world-line coordinate
$\gamma$, as obtained from the critical point of $\mathds{E}_{\rm IVP}$ with
$V(x)=\alpha x$.

Let us plot next in fig. 3, the results from our geometrized formalism as
physical trajectory, i.e. as ${\bf x}_{1,2}(t_{1,2})$ (red circles and blue
crosses). This allows us to compare the outcome to the solution one would
obtain by following the conventional approach in the literature (see e.g.
chapter 7.9 in goldstein1980classical ). There one considers time as
independent variable and simply adds a potential term to the free relativistic
action eq. 2 before deriving the corresponding Euler-Lagrange equation, which
for the linear potential reads
$d^{2}x/dt^{2}=-(\alpha)(1-(dx/dt)^{2})^{(3/2)}$. Using the LSODA algorithm of
Mathematica’s NDSolve command, we compute the solution of this equation of
motion and plot it as gray solid line. Excellent agreement with the solution
from our variational approach is observed, indicating that the geometrization
strategy indeed reproduces the solution of the physical problem at hand.

Note that the change in the velocity of the time coordinate manifests itself
here as a slightly denser time grid around the maximum of the trajectory.

Figure 3: The forward ${\bf x}_{1}(t_{1})$ (red circles) and backward ${\bf
x}_{2}(t_{2})$ (blue crosses) degrees of freedom obtained from the critical
point of $\mathds{E}^{\rm lin}_{\rm IVP}$ with $V(x)=\alpha x$. We discretize
with $N_{\gamma}=32$ points and the SBP21 operators. The trajectory from a
LSODA solver of the corresponding equation of motion is given as gray solid
line.

After this qualitative visual inspection, let us take a closer look at the
properties of the obtained solution. The first question we may ask is how well
quantitatively the solution follows the naively discretized geodesic equations
for time eq. 18 and space eq. 19 respectively. The continuum geodesic
equations for the system at hand read

$\displaystyle\frac{d}{d\gamma}\Big{(}g_{00}\frac{dt}{d\gamma}\Big{)}=\frac{d}{d\gamma}\Big{(}\big{(}1+2\alpha
x\big{)}\frac{dt}{d\gamma}\Big{)}=0,$ (70)
$\displaystyle\frac{d}{d\gamma}\Big{(}\frac{dx}{d\gamma}\Big{)}+\frac{1}{2}\frac{\partial
g_{00}}{\partial
x}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}=\frac{d^{2}x}{d\gamma^{2}}+\alpha\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}=0.$
(71)

When deriving these equations of motion from the continuum action functional
eq. 28 we have only used integration by parts. This motivates us to proceed,
considering them naively discretized by replacing the derivatives with SBP
finite difference operators

$\displaystyle\mathds{D}\big{(}(1+2\alpha{\bf x})\circ(\mathds{D}{\bf
t})\big{)}=\Delta{\bf G}^{t},$ (72) $\displaystyle\mathds{D}\mathds{D}{\bf
x}+\alpha(\mathds{D}{\bf t})\circ(\mathds{D}{\bf t})=\Delta{\bf G}^{x}.$ (73)

Here element-wise multiplication of entries of vector quantities is explicitly
denoted by the symbol $\circ$, which implements e.g. ${\bf x}_{1}\circ{\bf
x}_{2}=(x_{1}(0)x_{2}(0),x_{1}(\Delta\gamma)x_{2}(\Delta\gamma),\ldots,x_{1}(\Delta\gamma(N_{\gamma}-1))x_{2}(\Delta\gamma(N_{\gamma}-1)))^{\rm
T}$. Note that we have introduced on the right of the above equations two
quantities $\Delta{\bf G}^{x}$ and $\Delta{\bf G}^{t}$, which denote the
deviation from the value zero, to which the equations of motion evaluate in
the continuum. By inspecting $\Delta{\bf G}^{x}$ and $\Delta{\bf G}^{t}$ for
the trajectories ${\bf x}_{1,2}$ and ${\bf t}_{1,2}$ obtained from the
critical point of the discretized action functional $\mathds{E}^{\rm lin}_{\rm
IVP}$, we can obtain first quantitative insight into the performance of our
variational approach.

We plot the values of both quantities $\Delta{\bf G}^{x}$ and $\Delta{\bf
G}^{t}$ in the top panel of fig. 4. At first sight we find that deviations
from the naively discretized geodesic equations are minute, except for the two
last points. Note that the plot is given in logarithmic scale.

Since we use a minimizer in Mathematica with WorkingPrecision set to 40, the
values of $<10^{-30}$ reflect a true zero. It is apparent that both the
naively discretized geodesic equation for $x$ and $t$ are fulfilled down to
machine precision.

Figure 4: (top) Deviation of the spatial coordinate $\Delta{\bf G}^{x}$ (red
circles) and time coordinate $\Delta{\bf G}^{t}$ (blue crosses) from the
discretized geodesic equation, as obtained from the critical point of
$\mathds{E}^{\rm lin}_{\rm IVP}$ with $V(x)=\alpha x$ at $N_{\gamma}=32$
discretized with the SBP21 operator. (bottom) Deviation $\Delta{\bf E}$ of the
quantity ${\mathds{Q}}_{t}$ from its continuum value as given by the initial
conditions. Results for the SBP21 operator are given as red circles and those
for SBP42 as blue crosses. Note that the deviation $\Delta E$ in the interior
is exactly zero within machine precision.

Let us proceed to the central quantity of interest in this study $Q_{t}$,
defined in eq. 22, which in the continuum represents the conserved quantity
associated with the time-translation symmetry of the system. We again consider
its naively discretized form in the following

$\displaystyle{\bf\mathds{Q}_{t}}=(\mathds{D}{\bf t})\circ({\bf 1}+2\alpha{\bf
x}).$ (74)

With the discrete action functional $\mathds{E}^{\rm lin}_{\rm IVP}$ retaining
manifest invariance under shifts in the time coordinates ${\bf t}_{1,2}$ we
wish to investigate whether also the discretized $\mathds{Q}_{t}$ retains its
role as conserved Noether charge. To this end let us focus here on the
deviation $\Delta{\bf E}$ of $\mathds{Q}_{t}$ from its continuum value

$\displaystyle\Delta{\bf E}={\bf\mathds{Q}_{t}}-Q_{t}=(\mathds{D}{\bf
t})\circ({\bf 1}+2\alpha{\bf x})-\dot{t}_{i}(1+2\alpha x_{i}).$ (75)

Note that $\mathds{Q}_{t}$ takes on the continuum value by construction at the
first point in $\gamma$, as there it is defined by the initial conditions. The
values obtained for $\Delta{\bf E}$ from the critical point of
$\mathds{E}^{\rm lin}_{\rm IVP}$ using either the SBP21 (red circles) or SBP42
operator (blue crosses) are shown in the bottom panel of fig. 4. There are two
important observations to be made.

First, the discretized quantity $\mathds{Q}_{t}$ is exactly conserved in the
discrete setting in the interior of the simulated time domain and only at the
final point $\gamma_{f}$ it deviates from that constant. While the deviation
$\Delta{\bf E}(\gamma_{f})$ in case of the SBP21 operator is already smaller
than two permille, it reduces even further to a value of $10^{-6}$ when
deploying the SBP42 operator.

We have investigated various potential reasons for the slight difference at
the final point, such as a potential over-constraint from the connecting
conditions in eqs. 31 and 35, but we have not identified the source as of yet.
One venue to explore in the future is whether the exact enforcement of the
connecting conditions plays a role, which however requires the development of
a genuinely weak formulation of our approach without the use of Lagrange
multipliers. It is important to point out that, as we will show explicitly
below, the presence of this final differing point does not spoil the
convergence to the correct continuum limit.

Secondly, the value of $\mathds{Q}_{t}$ that remains conserved in the interior
agrees with the true continuum value, prescribed by the initial conditions,
within machine precision. This is a highly non-trivial result, as even in
energy preserving schemes, such as the leap-frog, the conserved quantities do
not necessarily agree with the continuum ones.

Figure 5: Visualization of the convergence of our novel variational approach
under the refinement of the $\gamma$ grid, based on the final point via
$\epsilon(\gamma_{f})$ (top panel) and globally via $\epsilon(L_{2})$ (bottom
panel). In each plot we show the results from the SBP21 (top two lines) and
SBP42 (bottom two lines) operators. Red circles and blue crosses denote the
absolute deviation in the spatial coordinate $\epsilon_{x}$ and temporal
coordinate $\epsilon_{t}$ between the continuum solution and the values
obtained from the the critical point of $\mathds{E}^{\rm lin}_{\rm IVP}$ with
$V(x)=\alpha x$. The corresponding deviations for the SBP42 operator are given
by the green circles and orange triangles respectively.

We surmise that it is the interplay of a manifest time-translation invariant
formulation of the action functional, together with the resulting dynamically
emerging time discretization, which achieves the conservation of the discrete
$\mathds{Q}_{t}$ at its continuum value in the interior of the simulation
domain.

The presence of two points that deviate from the naively discretized continuum
geodesic equations may appear troublesome. However as we show in fig. 5 these
points do not spoil the convergence to the correct continuum limit under grid
refinement.

In the top panel of fig. 5, we select the apparently most disadvantageous
points for our convergence study, i.e. we compare the deviation from the
continuum geodesic equations $\epsilon(\gamma_{f})_{x}=|{\bf
x}[N_{\gamma}]-x_{\rm true}(1)|$ and $\epsilon(\gamma_{f})_{t}=|{\bf
t}[N_{\gamma}]-t_{\rm true}(1)|$ at $\gamma_{f}$, exactly where the deviation
from the continuum result was maximal in the top panel of fig. 4. Grid
refinement is carried out and we provide the results for both the lowest order
SBP21 operator and the next higher SBP42 operator.

Even in this disadvantaged scenario, we find that under grid refinement, the
discrete solution approaches the true continuum values as expected from a
scheme that is second order in the interior. Taking the SBP21 results, the
best fit to $\epsilon_{x}$ reveals a scaling with $\Delta\gamma^{2.08}$, while
for $\epsilon_{t}$ an virtually identical $\Delta\gamma^{2.07}$ ensues. Going
over to the SBP42 results we find that the convergence is in line with
expectations for an SBP operator of 4th order in the interior with
$\epsilon_{x}$ exhibiting a scaling of $\Delta\gamma^{3.07}$ and
$\epsilon_{t}$ a somewhat better value of $\Delta\gamma^{3.48}$.

In the bottom panel of fig. 5 we instead investigate the global convergence of
our approach using the $L_{2}$ norm $\epsilon(L_{2})_{x}=\sqrt{({\bf x}-{\bf
x}_{\rm true})^{\rm T}.\mathds{H}.({\bf x}-{\bf x}_{\rm true})}$ and
$\epsilon(L_{2})_{t}=\sqrt{({\bf t}-{\bf t}_{\rm true})^{\rm
T}.\mathds{H}.({\bf t}-{\bf t}_{\rm true})}$, where ${\bf x}_{\rm true}$ and
${\bf t}_{\rm true}$ are taken from the numerical solution of the geodesic
equations, used for comparison in fig. 3. We find that similar convergence
rates ensue, where SBP21 shows scaling $\Delta\gamma^{\beta}$ with exponent
$\beta\geq 2$ and SBP42 shows scaling with exponent $\beta\geq 3$.

These convergence result agrees with the findings of our previous study
Rothkopf:2022zfb , where the standard action functional was discretized with
time as independent parameter.

###  4.2 Quartic potential

After considering the simplest possible non-trivial scenario with a linear
potential, we now turn to a system with a quartic potential and the following
continuum action functional

$\displaystyle{\cal E}^{\rm qrt}_{\rm IVP}=$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\frac{1}{2}\left\\{\Big{(}1+2\kappa
x_{1}^{4}(\gamma)\Big{)}\Big{(}\frac{dt_{1}}{d\gamma}\Big{)}^{2}-\Big{(}\frac{dx_{1}}{d\gamma}\Big{)}^{2}\right\\}$
$\displaystyle-$
$\displaystyle\int_{\gamma_{i}}^{\gamma_{f}}d\gamma\frac{1}{2}\left\\{\Big{(}1+2\alpha
x_{2}^{4}(\gamma)\Big{)}\Big{(}\frac{dt_{2}}{d\gamma}\Big{)}^{2}-\Big{(}\frac{dx_{2}}{d\gamma}\Big{)}^{2}\right\\}$
$\displaystyle+$
$\displaystyle\lambda_{1}\big{(}t(\gamma_{i})-t_{i}\big{)}+\lambda_{2}\big{(}\dot{t}_{1}(\gamma_{i})-\dot{t}_{i}\big{)}+\lambda_{3}\big{(}x_{1}(\gamma_{i})-x_{i}\big{)}+\lambda_{4}\big{(}\dot{x}_{1}(\gamma_{i})-\dot{x}_{i}\big{)}$
$\displaystyle+$
$\displaystyle\lambda_{5}\big{(}t_{1}(\gamma_{f})-t_{2}(\gamma_{2})\big{)}+\lambda_{6}\big{(}\dot{t}_{1}(\gamma_{f})-\dot{t}_{2}(\gamma_{2})\big{)}$
$\displaystyle+$
$\displaystyle\lambda_{7}\big{(}x_{1}(\gamma_{f})-x_{2}(\gamma_{2})\big{)}+\lambda_{8}\big{(}\dot{x}_{1}(\gamma_{f})-\dot{x}_{2}(\gamma_{2})\big{)}.$
(76)

Again we discretize along $N_{\gamma}=32$ in the world-line parameter
$\gamma$. Using $\kappa=1/2$ in the potential $V(x)=\kappa x^{4}$ leads to
dynamics that already in the small time regime considered here are distinctly
anharmonic.

Figure 6: The values of the (top) time coordinates ${\bf t}_{1}(\gamma_{i})$
(red circles) and ${\bf t}_{2}(\gamma_{i})$ (blue crosses) and that of the
(bottom) spatial coordinates ${\bf x}_{1}(\gamma_{i})$ (red circles) and ${\bf
x}_{2}(\gamma_{i})$ (blue crosses) along the world-line parameter $\gamma$, as
obtained from the critical point of $\mathds{E}^{\rm qrt}_{\rm IVP}$ with
$V(x)=\kappa x^{4}$, discretized with $N_{\gamma}=32$ points and the SBP21
operators. The solution of the corresponding geodesic equations via
Mathematica’s NDSolve is shown as light gray solid line.

As in the previous subsection we discretize the world-line of the particle
motion between $\gamma_{i}=0$ and $\gamma_{f}=1$, set the starting time to
$t_{i}=0$ and the starting position to $x_{i}=1$. For our choice of
$v_{i}=1/10$ we again decide on $\dot{t}=1$ and $\dot{x}=v_{i}$. The
discretized action functional thus reads

$\displaystyle\mathds{E}^{\rm qrt}_{\rm IVP}=$
$\displaystyle\frac{1}{2}\left\\{({\bar{\mathds{D}}^{\rm R}}_{t}{\bf
t}_{1})^{\rm T}\mathbb{d}\left[1+2\kappa{\bf
x}^{4}_{1}\right]{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{t}{\bf
t}_{1})-({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{1})^{\rm
T}{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{1})\right\\}$
$\displaystyle-$ $\displaystyle\frac{1}{2}\left\\{({\bar{\mathds{D}}^{\rm
R}}_{t}{\bf t}_{2})^{\rm T}\mathbb{d}\left[1+2\kappa{\bf
x}^{4}_{2}\right]{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{t}{\bf
t}_{2})-({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{2})^{\rm
T}{\bar{\mathds{H}}}({\bar{\mathds{D}}^{\rm R}}_{x}{\bf x}_{2})\right\\}$
$\displaystyle+$ $\displaystyle\lambda_{1}\big{(}{\bf
t}_{1}[1]-t_{i}\big{)}+\lambda_{2}\big{(}(\mathds{D}{\bf
t}_{1})[1]-\dot{t}_{i}\big{)}$ $\displaystyle+$
$\displaystyle\lambda_{3}\big{(}{\bf
x}_{1}[1]-x_{i}\big{)}+\lambda_{4}\big{(}(\mathds{D}{\bf
x}_{1})[1]-\dot{x}_{i}\big{)}$ $\displaystyle+$
$\displaystyle\lambda_{5}\big{(}{\bf t}_{1}[N_{\gamma}]-{\bf
t}_{2}[N_{\gamma}]\big{)}+\lambda_{6}\big{(}{\bf x}_{1}[N_{\gamma}]-{\bf
x}_{2}[N_{\gamma}]\big{)}$ $\displaystyle+$
$\displaystyle\lambda_{7}\big{(}(\mathds{D}{\bf
t}_{1})[N_{\gamma}]-(\mathds{D}{\bf
t}_{2})[N_{\gamma}]\big{)}+\lambda_{8}\big{(}(\mathds{D}{\bf
x}_{1})[N_{\gamma}]-(\mathds{D}{\bf x}_{2})[N_{\gamma}]\big{)}$ (77)

and taking the fourth power of the ${\bf x}_{1,2}$ vector is to be understood
in an element wise fashion.

While for the linear potential, the time geodesic appeared to depend almost
linearly on $\gamma$, we find that here a distinct curvature along $\gamma$
emerges, as shown in the top panel of fig. 6. We plot the values of ${\bf
t}_{1}(\gamma_{i})$ as red circles and ${\bf t}_{2}(\gamma_{i})$ as blue
crosses and show as gray solid line the solution of the corresponding geodesic
equation, obtained from the LSODA algorithm of Mathematica’s NDSolve command.
Again the physical limit of equal values ${\bf t}_{1}(\gamma)={\bf
t}_{2}(\gamma)$ is realized.

The values of the spatial coordinate ${\bf x}_{1}(\gamma_{i})$ and ${\bf
x}_{2}(\gamma_{i})$ as obtained from the critical point of $\mathds{E}^{\rm
qrt}_{\rm IVP}$ with $V(x)=\kappa x^{4}$ are plotted in the bottom panel of
fig. 6 with the direct numerical solution of the geodesic equation added as
gray solid line.

Figure 7: The change of time with respect to the world-line coordinate
$\gamma$, as obtained from the critical point of $\mathds{E}^{\rm qrt}_{\rm
IVP}$ with $V(x)=\kappa x^{4}$.

Note that even though we have provided an initial velocity of the time along
$\gamma$ again with value $\dot{t}_{i}=1$, the final time reached by the
simulation now lies at ${\bf t}[N_{\gamma}]=1.47$. Similarly one finds that
that a dynamical discretization in $t$ emerges, which, as shown in fig. 7,
varies from the initial values $\dot{t}_{i}=1$ to $(\mathds{D}{\bf
t})[N_{\gamma}]=2.06$. This behavior can be understood when realizing that the
trajectory $x(t)$ in the non-linear case shows a stronger curvature close to
$t=0$ than at later times. I.e. we find again that the automatically generated
non-trivial mesh (through automatic AMR) for the time coordinate adapts to the
dynamics, by exhibiting a finer spacing at initial times.

Let us take a look at the results from our geometrized formalism as physical
trajectory in fig. 8, i.e. plotted as ${\bf x}_{1,2}(t_{1,2})$ (red circles
and blue crosses). They are compared to the solution of the conventional
equation of motion, obtained from treating time as independent variable
$d^{2}x/dt^{2}=-(4\kappa x^{3})(1-(dx/dt)^{2})^{(3/2)}$, computed via the
LSODA algorithm of Mathematica’s NDSolve command (gray solid line) in the
range $t\in[0,1]$. We find that within this range the solution from our
geometrized discrete approach shows excellent agreement. Note that due to the
non-equidistant emergent time discretization, the physical trajectory $x(t)$,
shown in fig. 8 extends beyond the point $t=1$.

Figure 8: The forward ${\bf x}_{1}(t_{1})$ (red circles) and backward ${\bf
x}_{2}(t_{2})$ (blue crosses) degrees of freedom obtained from the critical
point of $\mathds{E}^{\rm qrt}_{\rm IVP}$ with $V(x)=\kappa x^{4}$. We
discretize with $N_{\gamma}=32$ points and the SBP21 operators. The trajectory
from a LSODA solver of the corresponding equation of motion is given as gray
solid line.

As for the linear potential, let us investigate quantitatively the properties
of the trajectories ${\bf t}(\gamma_{i})$ and ${\bf x}(\gamma_{i})$ by
inserting them into the naively discretized geodesic equations. For the
quartic potential, the continuum geodesic equations for the temporal and
spatial coordinate read

$\displaystyle\frac{d}{d\gamma}\Big{(}g_{00}\frac{dt}{d\gamma}\Big{)}=\frac{d}{d\gamma}\Big{(}\big{(}1+2\kappa
x^{4}\big{)}\frac{dt}{d\gamma}\Big{)}=0,$ (78)
$\displaystyle\frac{d}{d\gamma}\Big{(}\frac{dx}{d\gamma}\Big{)}+\frac{1}{2}\frac{\partial
g_{00}}{\partial
x}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}=\frac{d^{2}x}{d\gamma^{2}}+4\kappa
x^{3}\Big{(}\frac{dt}{d\gamma}\Big{)}^{2}=0.$ (79)

Naively discretizing these equations by replacing derivatives with SBP
operators leads to the following discrete geodesic equations

$\displaystyle\mathds{D}\big{(}(1+2\kappa{\bf x}^{4})\circ\mathds{D}{\bf
t}\big{)}=\Delta{\bf G}^{t},$ (80) $\displaystyle\mathds{D}\mathds{D}{\bf
x}+(4\kappa{\bf x}^{3})\circ(\mathds{D}{\bf t})\circ(\mathds{D}{\bf
t})=\Delta{\bf G}^{x}.$ (81)

where again taking a power of the ${\bf x}_{1,2}$ vector is to be understood
in an element wise fashion. To evaluate how well the solution obtained from
the critical point of $\mathds{E}^{\rm qrt}_{\rm IVP}$ fulfills the naive
discretized geodesic equations we have again introduced the quantities
$\Delta{\bf G}^{t}$ and $\Delta{\bf G}^{x}$ above.

As shown in fig. 9 also here in the highly non-linear scenario, we find that
the values of both ${\bf x}$ (red circles) and ${\bf t}$ (blue crosses) follow
the discretized geodesic equations excellently, except for the last two
points.

Figure 9: (top) Deviation of the spatial coordinate $\Delta{\bf G}^{x}$ (red
circles) and time coordinate $\Delta{\bf G}^{t}$ (blue crosses) from the
discretized geodesic equation, as obtained from the critical point of
$\mathds{E}^{\rm qrt}_{\rm IVP}$ with $V(x)=\kappa x^{4}$ at $N_{\gamma}=32$
discretized with the SBP21 operator. (bottom) Deviation $\Delta{\bf E}$ of the
quantity ${\mathds{Q}}_{t}$ from its continuum value as given by the initial
conditions. Results for the SBP21 operator are given as red circles and those
for SBP42 as blue crosses. Note that the deviation $\Delta{\bf E}$ in the
interior is exactly zero within machine precision.

The most important question however remains whether in the non-linear
discretized system, the continuum quantity $Q_{t}$ from eq. 22 also remains
conserved. Its naively discretized counterpart here reads

$\displaystyle{\bf\mathds{Q}_{t}}=(\mathds{D}{\bf t})\circ({\bf 1}+2\kappa{\bf
x}^{4}),$ (82)

and we define its deviation from the continuum result via the difference

$\displaystyle\Delta{\bf E}={\bf\mathds{Q}_{t}}-Q_{t}=(\mathds{D}{\bf
t})\circ({\bf 1}+2\kappa{\bf x}^{4})-\dot{t}_{i}(1+2\kappa x_{i}^{4}),$ (83)

which we plot in the bottom panel of fig. 9 using the SBP21 operator (red
circles) and the SBP42 operator (blue crosses).

We find also in the case of a non-linear potential that ${\bf\mathds{Q}_{t}}$
is preserved exactly in the interior of the simulation time domain. Up to
machine precision its values in the interior also take on the correct
continuum value. Similar to what we saw in the linear case, the last point
deviates from the continuum value. It is reassuring to see that the absolute
deviation at $\gamma_{f}$ reduces already by an order of magnitude when going
from a SBP21 to an SBP42 operator.

One may now ask whether the deviation of $\Delta{\bf E}$ from its continuum
value at $\gamma_{f}$ is in some way related to the fact that we use
$N_{\gamma}=32$ points to discretize the world-line parameter. The answer is
negative, as demonstrated in fig. 10.

Figure 10: Deviation of the quantity $\mathds{Q}_{t}$ from its continuum value
$Q_{t}$ in case of a fixed final time and grid refinement in $N_{\gamma}$.
Note that the magnitude of the deviation at the final point monotonously
decreases with grid refinement.

Three different datasets are shown in fig. 10, where for fixed $\dot{t}_{i}$
the grid spacing in $\gamma$ is changed. The green triangles denote the
results for $\Delta E$ when using $N_{\gamma}=16$, the red circles
$N_{\gamma}=32$ and the blue crosses $N_{\gamma}=64$. We have confirmed
explicitly that in all cases the values of $Q_{t}$ are preserved up to machine
precision in the interior of the simulated time domain. It is indeed only the
last point that shows a deviation and we see that the absolute magnitude of
the deviation reduces as the grid is refined.

Figure 11: (top) Deviation of the quantity $\mathds{Q}_{t}$ from its continuum
value $Q_{t}$ when moving the latest time to larger values by specifying
different initial $\dot{t}_{i}$. (bottom) The corresponding simulated
trajectories ${\bf x}(t)$.

For the next test, we instead increase $N_{\gamma}$ together with
$\dot{t}_{i}$ to let the simulation proceed to larger values of time $t$. In
the top panel of fig. 11 we plot the deviation of $Q_{t}$ from its continuum
value for three choices $\dot{t}_{i}=1,N_{\gamma}=16$ (green triangles),
$\dot{t}_{i}=4,N_{\gamma}=32$ (red circles) and $\dot{t}_{i}=8,N_{\gamma}=64$
(blue crosses). As seen before in the interior of the simulated time domain,
the values of $\mathds{Q}_{t}$ remain exactly preserved and only the last
point deviates. We find that the magnitude of the deviation in the last point
changes only marginally with the length of the simulated trajectory. For
completeness the corresponding trajectories $x(t)$ are plotted in the bottom
panel of fig. 11. Again let us emphasize that, as we will show below, the
presence of this single deviating point does not spoil the convergence to the
correct solution under grid refinement.

The exact conservation of the quantity ${\mathds{Q}}_{t}$ in the interior is
remarkable, as e.g. the trajectory in the bottom panel of fig. 11 for
$\dot{t}_{i}=8,N_{\gamma}=64$ shows sizable discretization artifacts (which
disappear under grid refinement). We believe that it is due to the manifest
time-translation invariance of the underlying action functional that the
combined dynamics of $x(\gamma)$ and $t(\gamma)$, including the automatically
generated non-equidistant time mesh, achieve conservation of the continuum
quantity.

The fact that the solutions we obtain fulfill the naively discretized geodesic
equations and provide exact conservation of the continuum conserved charge in
the interior of the simulated domain (see fig. 9) bodes well for establishing
its stability. Since in the IVP setting $t(\gamma_{f})$ is not given but
emerges dynamically we cannot directly apply eq. 26 as proof of stability.
However, as long as we can assume that the simulated time range (given a
certain $\dot{t}(\gamma_{i})$ is finite, the linear bound of eq. 26 on the
norm ${\cal H}_{\rm BVP}$ holds in the discrete setting. In turn we deduce
that the solution cannot exhibit stronger than linear rise of the derivatives
of either $t(\gamma)$ or $x(\gamma)$, implying stability of the approach.

Let us now quantify the convergence properties of our variational approach
using the results from the lowest order SBP21 operator and those coming from
the SBP42 operator in fig. 12.

Figure 12: Visualization of the convergence of our novel variational approach
under the refinement of the $\gamma$ grid, based on the final point via
$\epsilon(\gamma_{f})$ (top panel) and globally via $\epsilon(L_{2})$ (bottom
panel). In each plot we show the results from the SBP21 (top two lines) and
SBP42 (bottom two lines) operators. Red circles and blue crosses denote the
absolute deviation in the spatial coordinate $\epsilon_{x}$ and temporal
coordinate $\epsilon_{t}$ between the continuum solution and the values
obtained from the the critical point of $\mathds{E}^{\rm qrt}_{\rm IVP}$ with
$V(x)=\kappa x^{4}$. The corresponding deviations for the SBP42 operator are
given by the green circles and orange triangles respectively.

As in the linear potential case, in the top panel of fig. 12, we select the
most disadvantageous points for our convergence study, i.e. we compare the
deviation from the continuum geodesic equations
$\epsilon(\gamma_{f})_{x}=|{\bf x}[N_{\gamma}]-x_{\rm true}(1)|$ and
$\epsilon(\gamma_{f})_{t}=|{\bf t}[N_{\gamma}]-t_{\rm true}(1)|$ at
$\gamma_{f}$, exactly where the deviation from the continuum result was
maximal in the top panel of fig. 9. Also in the non-linear scenario we find
that under refinement of the $\gamma$ grid, the discrete solution monotonously
approaches the true continuum values.

Taking the SBP21 results, the best fit to $\epsilon(\gamma_{f})_{x}$ reveals a
scaling with $\Delta\gamma^{2.08}$, while for $\epsilon(\gamma_{f})_{t}$ an
virtually identical $\Delta\gamma^{2.06}$ is obtained.

For SBP42, we find that the convergence is slightly worse than in the linear
potential case. As seen in the green circles plotted in fig. 12, the
asymptotic convergence regime is reached for $32<N_{\gamma}<64$. Once we are
in that regime, we find that $\epsilon(\gamma_{f})_{x}$ exhibits a scaling of
$\Delta\gamma^{2.84}$, close to the expected value of three. On the other hand
$\epsilon(\gamma_{f})_{t}$ shows a consistent performance with a scaling of
$\Delta\gamma^{3.13}$ already at $N_{\gamma}=32$.

Let us now investigate the global convergence in the bottom panel of fig. 12
using the $L_{2}$ norm $\epsilon(L_{2})_{x}=\sqrt{({\bf x}-{\bf x}_{\rm
true})^{\rm T}.\mathds{H}.({\bf x}-{\bf x}_{\rm true})}$ and correspondingly
$\epsilon(L_{2})_{t}=\sqrt{({\bf t}-{\bf t}_{\rm true})^{\rm
T}.\mathds{H}.({\bf t}-{\bf t}_{\rm true})}$, where ${\bf x}_{\rm true}$ and
${\bf t}_{\rm true}$ are taken from the numerical solution of the geodesic
equations, used for comparison in fig. 6.

Reassuringly we find that the global convergence properties of our approach
are better than indicated by those of the most disadvantaged point in the top
panel of fig. 12. Indeed we find that for the SBP41 operators, the global
scaling regime is reached already at $N_{\gamma}=32$, similarly to the SBP21
case. In addition, the global convergence rate $\Delta\gamma^{\beta}$ for
SBP42 operators lies consistently above $\beta\geq 3$ for both the $x$ and $t$
degrees of freedom.

Again, these convergence result are in good agreement with those of our
previous study Rothkopf:2022zfb , where the standard action functional was
discretized with time as independent parameter.

##  5 Summary and Outlook

In this study we have put forward a novel geometric variational approach for
solving a large class of initial value problems, associated with the dynamics
of point particles evolving under a generic $x$ dependent potential $V(x)$.
Taking inspiration from the general theory of relativity, we consider both
time and spatial coordinates of the point particle as dependent variables of a
world-line parameter $\gamma$. We select a continuum action functional, which
in the non-relativistic limit reduces to the standard action of point
mechanics and whose critical point encodes a set of geodesic equations for
$x(\gamma)$ and $t(\gamma)$. After doubling the degrees of freedom $t_{1,2}$
and $x_{1,2}$ we can relate the critical point of the corresponding doubled
d.o.f. action with the classical trajectory. Using the concept of Killing
vectors we identify conserved quantities, e.g. related to the continuum time
translation invariance of the action.

Deploying the regularized SBP operators originally introduced in
Rothkopf:2022zfb , we discretize the continuum action and add Lagrange
multipliers to enforce the initial and connecting conditions between the
doubled ${\bf t}_{1,2}$ and ${\bf x}_{1,2}$. The main novelty of our approach
is that the discretized action retains the continuum symmetries, in particular
the invariance under time translations. Exactly mimicking integration by part
through the use of SBP finite difference operators entails that the derivation
of the conserved charges associated with the Killing vectors of the system is
also exactly mimicked in the discrete setting. I.e. the continuum conserved
quantities $Q_{K}$ retain their role even after discretization.

The numerical results we obtain for both a linear and highly non-linear
potential show that a discretization of time $t$ now indeed emerges
dynamically, adapting to the behavior of the spatial coordinate $x$. This is a
concrete realization of an automatically generated non-equidistant mesh for
the time coordinate, guided by our action functional with manifest continuum
translation symmetry, i.e. an automatic AMR procedure. We have shown that
except for the last two points along the discrete $\gamma$, the solution we
obtain follows the naively discretized geodesic equations excellently.

Even more importantly, the naively discretized counterpart $\mathds{Q}_{t}$ of
the continuum conserved quantity $Q_{t}$ remains exactly preserved in the
interior of the simulated time domain, where it even retains its continuum
value exactly within machine precision. A small deviation from the values in
the interior for $\mathds{Q}_{t}$ is observed at the last step $\gamma_{f}$.
This deviation however decreases both under grid refinement, as well as when
increasing the order of the SBP operator.

Point-wise, as well as global scaling analyses under grid refinement show that
even in the presence of two points deviating from the naively discretized
geodesic equations at the last two $\gamma$ steps, the solution monotonously
improves and manages to approach the true solution. When deploying the SBP21
operator, we achieve consistent scaling in $\Delta\gamma^{\beta}$ with
$\beta\gtrsim 2$ for both the linear and non-linear potential. For SBP42 in
case of a linear potential the dependence on the grid spacing follows the
expected power law $\Delta\gamma^{\beta}$ with $\beta\gtrsim 3$ for all values
of $N_{\gamma}$ we inspected. For the non-linear potential, the scaling regime
for point-wise convergence at the last point $\gamma_{f}$ is reached with
SBP42 for $32<N_{\gamma}<64$ with a slightly worse scaling of
$2.84\leq\beta\leq 3.13$. Global convergence on the other hand shows
consistent scaling at all $N_{\gamma}$ we considered, with exponents
$\beta\geq 3$, in agreement with the findings in our previous paper
Rothkopf:2022zfb , where the standard action functional was discretized with
time as independent variable.

This study presents a proof of principle that initial value problems can be
discretized, while retaining continuum symmetries. Three future directions
will be explored: we may ask how we can capture systems of ordinary
differential equations that e.g. contain a term that is proportional to a
first derivative in $x$ with respect to time? To this end we must exploit the
versatility of the doubled d.o.f. approach more thoroughly. Furthermore we
will explore how the reparametrization invariant formulation can be applied to
partial differential equations in higher dimensions, taking insight from how
the non-relativistic action emerges from our relativistic starting point in
eq. 8. In addition, to better understand the origin of the single deviating
value in the otherwise exactly preserved $\mathds{Q}_{t}$, we will develop a
genuinely weak formulation of our approach, devoid of Langrange multipliers
for enforcing initial and connecting conditions.

We believe that the quest for retention of defining continuum properties in
discretized systems is both conceptually and practically valuable. Not only
does the preservation of symmetries place powerful physical constraints on the
solution but in addition offers a mechanism for the automatic generation of
optimal discrete spacetime grids to ensure conservation of the Noether charges
associated with these symmetries. We hope that this study provides the
community with a novel impulse in this direction.

## Acknowledgements

A. R. thanks Will Horowitz for inspiring and insightful discussions and Alex
Nielsen for valuable insight on the general theory of relativity. A. R. gladly
acknowledges support by the Research Council of Norway under the FRIPRO Young
Research Talent grant 286883. J. N. was supported by the Swedish Research
Council grant nr. 2021-05484. The study has benefited from computing resources
provided by UNINETT Sigma2 - the National Infrastructure for High Performance
Computing and Data Storage in Norway under project NN9578K-QCDrtX ”Real-time
dynamics of nuclear matter under extreme conditions”

##  Appendix A Regularized SBP operators in affine coordinates

We here briefly review the idea and some technical aspects of constructing
null-space consistent regularized SBP operators using affine coordinates,
developed in our study Rothkopf:2022zfb .

The goal of regularizing conventional SBP operators $\mathds{D}$, such as
those defined e.g. in eq. 38 and eq. 52, lies in removing their unphysical
zero modes. These may appear as highly oscillatory eigenfunctions to
$\mathds{D}^{\rm T}$ with zero eigenvalue. To this end we take inspiration
from regularization techniques developed for partial differential equations.
There the concept of null-space consistent SBP operators has been discussed in
detail (see e.g. svard2019convergence ; linders_properties_2020 ;
svard_convergence_2021 ; ranocha2021new ).

For a differential equation, the boundary conditions may be enforced in the
weak sense by adding a simultaneous approximation penalty term (SAT)
carpenter1994time , which can be partially absorbed into the finite difference
operator, lifting its zero modes. Take for example a simple discretized first
order differential equation

$\displaystyle\mathds{D}{\bf u}=\lambda{\bf
u}+\sigma_{0}\mathds{H}^{-1}\mathds{E}_{0}\big{(}{\bf u}-{\bf g}\big{)},$ (84)

where the SAT penalty term has been added to the right-hand side. It features
the matrix $\mathds{E}_{0}={\rm diag}[1,0,\ldots,0]$ that makes reference only
to the first entry in the discretized functions ${\bf u}$ and ${\bf g}$, the
latter of which contains the initial value in its first entry ${\bf
g}=(u_{0},0,\ldots,0)$. The SAT term also contains $\mathds{H}^{-1}$, i.e.
${\Delta t}^{-1}$, which increases the strength of the penalty as ${\Delta
t}\to 0$. The parameter $\sigma_{0}$ in the SBP-SAT approach is tuned to
satisfy stability properties and its optimal value is found to be
$\sigma_{0}=-1$ (see e.g. ref. aalund2016provably ; ALUND2019209 ;
glaubitz2023summation ), a choice we adopt in the following. In the
differential equation context one conventionally absorbs the term proportional
to ${\bf u}$ into a new
$\tilde{\mathds{D}}=\mathds{D}-\sigma_{0}\mathds{H}^{-1}E_{0}$. This new
operator is devoid of zero modes ruggiu_eigenvalue_2020  and may be inverted
to obtain the solution ${\bf u}$.

In the context of an action functional, such as eq. 67, we do not have an
equal sign around which we can move the SAT term. Instead we must incorporate
the whole of the penalty term directly in a modified SBP operator. Since the
penalty term in our example eq. 84 contains both a contribution that is
proportional to the function ${\bf u}$ and a constant shift ${\bf g}$ it
amounts to an affine transformation on ${\bf u}$, which can be captured
efficiently using affine coordinates. To this end let us write $\bar{A}[{\bf
b}]\bar{\bf x}=A{\bf x}+{\bf b}$, where $\bar{A}[{\bf b}]$ refers to a matrix
$A$ extended by an additional row and column with the value $1$ placed in the
lower right corner. The new column available in $\bar{A}[{\bf b}]$ is
populated with the entries of ${\bf b}$. The vector $\bar{\bf x}$ is nothing
but ${\bf x}$ extended by one more entry with value unity. We will use this
construction principle to define a regularized $\bar{\mathds{D}}$ from our
conventional SBP operator $\mathds{D}$.

Since we have both ${\bf x}$ and ${\bf t}$ as independent degrees of freedom
each with independent initial conditions $x_{i}$ and $t_{i}$, we must define
different shifts ${\bf b}^{x}$ and ${\bf b}^{t}$ respectively and thus end up
with two different regularized SBP operators ${\bar{\mathds{D}}}_{t}$ and
${\bar{\mathds{D}}}_{x}$. The shift terms are nothing but the constant part of
the corresponding SAT term, absorbed into the SBP operator

$\displaystyle{\bf b}^{x}=\sigma_{0}\mathds{H}^{-1}E_{0}{\bf g}^{x},\quad{\bf
b}^{t}=\sigma_{0}\mathds{H}^{-1}E_{0}{\bf g}^{t}.$ (85)

Here ${\bf g}^{x}={\rm diag}[x_{i},0,\cdots,0]$ and ${\bf g}^{t}={\rm
diag}[t_{i},0,\cdots,0]$ encode the initial values for $x$ and $t$
respectively. As mentioned before, we choose the parameter $\sigma_{0}=-1$,
whenever a penalty term is incorporated in ${\bar{\mathds{D}}}$, motivated by
the fact that in the conventional treatment of IVPs using the SBP-SAT
approach, this value leads to a minimal discretization error (see e.g. ref.
aalund2016provably ; ALUND2019209 ; glaubitz2023summation ). The resulting
regularized SBP operators to be deployed on ${\bf t}_{1,2}$ or ${\bf
x}_{1,2}$, are given explicitly in eq. 59 and eq. 66 respectively.

Consistent with the affine coordinates used in the newly defined
${\bar{\mathds{D}}}_{t}$ and ${\bar{\mathds{D}}}_{x}$, we also amend the
discretized trajectories ${\bf t}_{1,2}$ and ${\bf x}_{1,2}$ by one more entry
that is given the value one.

In order to compute inner products in the space of discretized functions, we
also have to modify the quadrature matrix $\mathds{H}\to\bar{\mathds{H}}$ by
amending it by one row and column filled with zeros. We do not include the
value one in the lower right corner in order to correctly account for the fact
that the vectors appearing as arguments to the inner product contain an
auxiliary final entry, which does not contribute to the value of the inner
product and only facilitates the efficient implementation of shift operations.
For more details on the affine coordinate regularization technique see
Rothkopf:2022zfb .

## Competing interests

The authors declare that they have no competing interests.

## Author’s contributions

  * •

A. Rothkopf: formulation of the geometric variational approach, literature
review, numerical experiments, writing, editing

  * •

J. Nordstöm: guidance on the formulation and implementation of SBP based
discretization schemes, literature review, editing

## References

  * (1) Goldstein, H., Poole, C.P., Safko, J.L.: Classical Mechanics. Addison Wesley (2002)
  * (2) Arnold, V.I., Vogtmann, K., Weinstein, A.: Mathematical Methods of Classical Mechanics. Graduate Texts in Mathematics. Springer (2013)
  * (3) Coleman, S.: Aspects of Symmetry: Selected Erice Lectures. Cambridge University Press, Cambridge, U.K. (1985). doi:10.1017/CBO9780511565045
  * (4) Noether, E.: Invariant variation problems. Transport theory and statistical physics 1(3), 186–207 (1971)
  * (5) Landau, L.D., Lifshitz, E.M.: The Classical Theory of Fields: Volume 2. Course of theoretical physics. Elsevier Science (2000)
  * (6) Yanagihara, R., Iritani, T., Kitazawa, M., Asakawa, M., Hatsuda, T.: Distribution of Stress Tensor around Static Quark–Anti-Quark from Yang-Mills Gradient Flow. Phys. Lett. B 789, 210–214 (2019). doi:10.1016/j.physletb.2018.09.067. 1803.05656
  * (7) Cockburn, B., Karniadakis, G.E., Shu, C.-W.: Discontinuous Galerkin Methods: Theory, Computation and Applications vol. 11. Springer (2012)
  * (8) Svärd, M., Nordström, J.: Review of summation-by-parts schemes for initial–boundary-value problems. Journal of Computational Physics 268, 17–38 (2014)
  * (9) Fernández, D.C.D.R., Hicken, J.E., Zingg, D.W.: Review of summation-by-parts operators with simultaneous approximation terms for the numerical solution of partial differential equations. Computers & Fluids 95, 171–196 (2014)
  * (10) Lundquist, T., Nordström, J.: The SBP-SAT technique for initial value problems. Journal of Computational Physics 270, 86–104 (2014)
  * (11) Nordström, J., Lundquist, T.: Summation-by-parts in time. Journal of Computational Physics 251, 487–499 (2013)
  * (12) Nordström, J., Lundquist, T.: Summation-by-parts in time: the second derivative. SIAM Journal on Scientific Computing 38(3), 1561–1586 (2016)
  * (13) Johnson, R.C.: Angular momentum on a lattice. Physics Letters B 114(2), 147–151 (1982). doi:10.1016/0370-2693(82)90134-4
  * (14) Regan, H.M.: Von neumann stability analysis of symplectic integrators applied to hamiltonian pdes. Journal of Computational Mathematics 20(6), 611–618 (2002)
  * (15) Nordström, J.: Nonlinear Boundary Conditions for Initial Boundary Value Problems with Applications in Computational Fluid Dynamics. 2306.01297
  * (16) Verlet, L.: Computer ”Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Physical Review 159(1), 98–103 (1967). doi:10.1103/PhysRev.159.98. Publisher: American Physical Society
  * (17) Dirac, P.a.M.: Generalized Hamiltonian Dynamics. Canadian Journal of Mathematics 2, 129–148 (1950). doi:10.4153/CJM-1950-012-1. Publisher: Cambridge University Press
  * (18) Anerot, B., Cresson, J., Hariz Belgacem, K., Pierret, F.: Noether’s-type theorems on time scales. Journal of Mathematical Physics 61(11), 113502 (2020). doi:10.1063/1.5140201. Number: 11
  * (19) Stephani, H.: Relativity: An Introduction to Special and General Relativity. Cambridge University Press (2004)
  * (20) Rothkopf, A., Nordström, J.: A new variational discretization technique for initial value problems bypassing governing equations. J. Comput. Phys. 477, 111942 (2023). doi:10.1016/j.jcp.2023.111942. 2205.14028
  * (21) Berger, M.J., Oliger, J.: Adaptive mesh refinement for hyperbolic partial differential equations. Journal of computational Physics 53(3), 484–512 (1984)
  * (22) Löhner, R.: An adaptive finite element scheme for transient problems in CFD. Computer methods in applied mechanics and engineering 61(3), 323–338 (1987)
  * (23) Berger, M.J., Colella, P.: Local adaptive mesh refinement for shock hydrodynamics. Journal of computational Physics 82(1), 64–84 (1989)
  * (24) Persson, P.-O., Peraire, J.: Sub-cell shock capturing for discontinuous Galerkin methods. In: 44th AIAA Aerospace Sciences Meeting and Exhibit, p. 112 (2006)
  * (25) Nemec, M., Aftosmis, M., Wintzer, M.: Adjoint-based adaptive mesh refinement for complex geometries. In: 46th AIAA Aerospace Sciences Meeting and Exhibit, p. 725 (2008)
  * (26) Offermans, N., Massaro, D., Peplinski, A., Schlatter, P.: Error-driven adaptive mesh refinement for unsteady turbulent flows in spectral-element simulations. Computers & Fluids 251, 105736 (2023)
  * (27) Mavriplis, C.: Adaptive mesh strategies for the spectral element method. Computer methods in applied mechanics and engineering 116(1-4), 77–86 (1994)
  * (28) Henderson, R.D.: Adaptive spectral element methods for turbulence and transition. In: High-order Methods for Computational Physics, pp. 225–324. Springer (1999)
  * (29) Kompenhans, M., Rubio, G., Ferrer, E., Valero, E.: Adaptation strategies for high order discontinuous Galerkin methods based on tau-estimation. Journal of Computational Physics 306, 216–236 (2016)
  * (30) Galley, C.R.: Classical Mechanics of Nonconservative Systems. Physical Review Letters 110(17), 174301 (2013). doi:10.1103/PhysRevLett.110.174301. Publisher: American Physical Society
  * (31) Jost, J., Li-Jost, X.: Calculus of Variations. Cambridge Studies in Advanced Mathematics. Cambridge University Press (1998)
  * (32) Carroll, S.M.: Spacetime and Geometry. Cambridge University Press (2019)
  * (33) Carlip, S.: General Relativity: A Concise Introduction. OUP Oxford (2019)
  * (34) Rizzuti, B.F., Júnior, G.F.V., Resende, M.A.: To square root the Lagrangian or not: an underlying geometrical analysis on classical and relativistic mechanical models. arXiv. arXiv:1905.01177 [math-ph, physics:physics] (2019). doi:10.48550/arXiv.1905.01177. http://arxiv.org/abs/1905.01177
  * (35) Berges, J., Gasenzer, T.: Quantum versus classical statistical dynamics of an ultracold Bose gas. Phys. Rev. A 76, 033604 (2007). doi:10.1103/PhysRevA.76.033604. cond-mat/0703163
  * (36) Carpenter, M.H., Gottlieb, D., Abarbanel, S.: Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodology and application to high-order compact schemes. Journal of Computational Physics 111(2), 220–236 (1994)
  * (37) Rothkopf, A.: Mathematica 12 implementation of a symmetry and Noether charge preserving IVP discretization technique (2023). doi:10.5281/zenodo.8129657
  * (38) Svärd, M., Nordström, J.: On the convergence rates of energy-stable finite-difference schemes. Journal of Computational Physics 397, 108819 (2019)
  * (39) Linders, V., Nordström, J., Frankel, S.H.: Properties of Runge-Kutta-Summation-By-Parts methods. Journal of Computational Physics 419, 109684 (2020). doi:10.1016/j.jcp.2020.109684
  * (40) Svärd, M., Nordström, J.: Convergence of energy stable finite-difference schemes with interfaces. Journal of Computational Physics 429, 110020 (2021). doi:10.1016/j.jcp.2020.110020
  * (41) Ranocha, H., Nordström, J.: A new class of a stable summation by parts time integration schemes with strong initial conditions. Journal of Scientific Computing 87(1), 1–25 (2021)
  * (42) Ålund, O., Nordström, J.: A provably stable, non-iterative domain decomposition technique for the advection-diffusion equation. Linköping University Electronic Press (2016)
  * (43) Ålund, O., Nordström, J.: Encapsulated high order difference operators on curvilinear non-conforming grids. Journal of Computational Physics 385, 209–224 (2019). doi:10.1016/j.jcp.2019.02.007
  * (44) Glaubitz, J., Nordström, J., Öffner, P.: Summation-by-parts operators for general function spaces. SIAM Journal on Numerical Analysis 61(2), 733–754 (2023)
  * (45) Ruggiu, A.A., Nordström, J.: Eigenvalue Analysis for Summation-by-Parts Finite Difference Time Discretizations. SIAM Journal on Numerical Analysis 58(2), 907–928 (2020). doi:10.1137/19M1256294. Publisher: Society for Industrial and Applied Mathematics