Exact symmetry conservation and automatic mesh refinement in discrete initial boundary value problems

]ARA. Rothkopf

]WAHW. A. Horowitz

]JNJ. Nordström

Faculty of Science and Technology,
University of Stavanger,

Department of Physics,
University of Cape Town,
South Africa

Department of Mathematics,
Linköping University,
SE-581 83,

Department of Mathematics and Applied Mathematics,
University of Johannesburg,
P.O. Box 524, Auckland Park 2006,
South Africa

We present a novel solution procedure for initial boundary value problems. The procedure is based on an action principle, in which coordinate maps are included as dynamical degrees of freedom. This reparametrization invariant action is formulated in an abstract parameter space and an energy density scale associated with the space-time coordinates separates the dynamics of the coordinate maps and of the propagating fields. Treating coordinates as dependent, i.e. dynamical quantities, offers the opportunity to discretize the action while retaining all space-time symmetries and also provides the basis for automatic adaptive mesh refinement (AMR). The presence of unbroken space-time symmetries after discretization also ensures that the associated continuum Noether charges remain exactly conserved. The presence of coordinate maps in addition provides new freedom in the choice of boundary conditions. An explicit numerical example for wave propagation in $1+1$ dimensions is provided, using recently developed regularized summation-by-parts finite difference operators.

Initial Boundary Value Problems, Space-Time Symmetry, Automatic Mesh Refinement, Noether Charge

§ MOTIVATION

§.§ Executive Summary

In this study we set out to address three central challenges in the treatment of discretized initial boundary value problems (IBVPs):

* the breaking of space-time symmetries related to the finite grid discretization and the associated lack of conservation of continuum Noether charges,

* the need to construct appropriate meshes to accurately resolve the simulated dynamics, and

* the need for more flexible and less costly implementation  of (non-reflecting) boundary conditions.

We will work directly on the level of the action of the system, from which governing equations can be derived. Our novel action is formulated in a set of abstract temporal and spatial parameters $(\tau,\vec{\sigma})$. While in the usual IBVP treatment conventional space-time coordinates $(t,\vec{x})$ represent the independent parameters of the theory, here space-time coordinates are introduced as dynamical mappings $(t(\tau,\vec{\sigma}),\vec{x}(\tau,\vec{\sigma}))$, which themselves depend on abstract parameters. When discretizing the abstract parameters, the coordinate maps remain continuous and our action retains its continuum space-time symmetries after discretization. As a result, this novel discretization technique preserves the continuum Noether charges exactly. Concurrently we show that the presence of dynamical coordinate maps leads to new contributions in boundary terms, offering novel freedom to implement boundary conditions. In addition, the dynamical interplay of the coordinate maps and the fields allows the coordinates to adapt and to form a non-trivial space-time mesh guided by the symmetries of the system, realizing a form of automatic adaptive mesh refinement.

A sketch of the differences between the conventional approach which discretizes space-time coordinates directly and our novel approach, which discretizes the underlying abstract parameters is provided in <ref>.

(Top) Sketch of the conventional approach to IBVPs: space-time coordinates are designated as independent variables and the field $\phi(t,x)$ propagates on the background of this space-time scaffold. When discretizing space-time as a hypercubic grid, the field $\phi(t_i,x_j)$ is resolved in a similarly regular fashion. Since discretized space-time is unable to accommodate infinitesimal symmetry transformations, the continuum Noether charge is not preserved. (Bottom) Sketch of our novel approach to IBVPs: A set of abstract parameters $(\tau,\sigma)$ is designated as independent parameters. Both the field $\phi(\tau,\sigma)$ as well as dynamical coordinate maps $t(\tau,\sigma)$ and $x(\tau,\sigma)$ evolve on the background of the $(\tau,\sigma)$ parameters. The evolution of the coordinate maps can be highly non-linear depending on the field dynamics. Discretizing the $(\tau,\sigma)$ parameters leaves the values of the coordinate maps continuous. In turn when expressing the physical field solution in terms of space-time coordinates $\phi(t(\tau_i,\sigma_j),x(\tau_i,\sigma_j))$, one finds in general that a non equidistant space-time grid emerges, which automatically adapts in resolution to the dynamics of the field. Since the discrete action retains its continuum symmetries the continuum Noether charge remains exactly conserved.

As a preview to the reader, we list the main achievements of our novel approach to IBVPs in the inset box on the next page. It showcases concrete examples from a numerical investigation of wave propagation in $(1+1)$d dimensions, treated via our novel approach, as presented in detail in <ref>.

[width=1]Executive Summary - with results from (1+1)d wave propagation (c.f. <ref>)

* Space-time coordinates are promoted to dynamical degrees of freedom, which propagate together with the fields in an abstract parameter space.

0.75(left) Propagating wave packages reflecting from a fixed spatial boundary. (center) Evolution of the dynamical time mapping. (right) The corresponding spatial coordinate mapping. (c.f. <ref> in <ref>)

* Interplay between fields and coordinate maps leads to a coarser or finer space-time resolution, depending on where relevant changes occur in the field configuration. This constitutes a form of dynamic resolution of the space-time coordinates which realizes automatic adaptive mesh refinement. This mechanism is independent of the specific discretization scheme and whether the IBVP is solved on the level of the action or governing equations.

0.75Non-trivial evolution of the time-mapping visualized via its temporal derivative. Larger values indicate a coarser time grid, smaller values a finer time grid. (c.f. <ref> in <ref>)

* Discretization in the abstract parameters instead of in space-time allows us to keep the coordinate maps continuous and thus the discretized action retains its space-time symmetries.

* If a discretization scheme is used that exactly mimics integration by parts, the presence of continuum space-time symmetries in the discretized action allows us to establish a discrete Noether theorem and the associated Noether charges are exactly preserved.

0.75The exactly preserved Noether charge associate with time translation symmetry. (c.f. <ref> in <ref>)

* With both fields and coordinates entering as dynamical degrees of freedom, the boundary terms we encounter contain novel contributions, which provide more flexibility in the construction of boundary conditions.

§.§ Introduction

Progress in both fundamental and applied science relies to a large extent on the solution of initial boundary value problems (IBVPs). Wave propagation is a prominent and important example. In the form of electromagnetic waves (for a textbook see e.g. [1]) accurate numerical simulations of wave propagation are a vital ingredient in wireless telecommunication and stealth defence technology. Wave propagation is also indispensable for understanding the minute ripples in space-time emanating from the merger of two black holes in the form of gravitational waves [2]. In the context of mechanical waves, the propagation of sound in various media and environments has received great attention (see e.g. [3]). The behavior of fluids is yet another area where IBVPs take center stage, be it non-relativistic flow (for a textbook see e.g. [4]) over airfoils and through intricate tubing or the relativistic flow (see e.g. [5]) of the almost perfect fluid of hot nuclear matter created in relativistic heavy-ion collisions.

In fact, all of the fundamental forces: gravity, electromagnetism, weak- and strong nuclear force are described in the language of a field theory [6, 7, 8, 9]. Their classical dynamics are conventionally formulated in the form of IBVPs, i.e. a set of coupled partial differential equations supplied with initial and boundary data related to a specific experimental or observational scenario. From a theory point of view, the fundamental mathematical object to summarize the classical mechanics of these fields is their action [10], defined as the time integral over the Lagrange functional. Dirac [11] and Feynman's [12] vital contribution to quantum field theory showed that also in the microscopic realm, the dynamics of fields is governed by this classical action. The presence  of quantum fluctuations leads to deviations from the classical configuration, but it is still the classical action that steers those quantum excursions.

The treatment of IBVPs on the level of their governing equations and through their action is intimately related. Take point mechanics as an example (see [10] for an in-depth introduction). We may understand the motion of particles either by using Newton's second law to analyze forces or we may gain insight by analyzing the energy budget of the system at hand. Both approaches must and do in the end lead to the same results. Governing equations are often motivated by considering the processes in a system, which lead to a change of its state, similar to applying Newton's second law in mechanics. To construct an action on the other hand, requires us to inspect the imbalance between the kinetic and potential energy contributions to our system. It is Hamilton's variational principle that relates the two approaches [10]. Hamilton's principle states that the classical trajectory or field configuration corresponds to a critical point of the classical action. By performing variational analysis (for a mathematical perspective see [13]) on the action, one can show that under certain boundary conditions, the critical point of the action can also be obtained by solving a set of differential equations, the Euler-Lagrange equations. I.e. we may equivalently either solve the Euler-Lagrange equations or find the critical point of the system action. For multi-component systems, it turns out that it often requires less effort to formulate the action than to construct the corresponding governing equations.

Working with the action has a clear theoretical advantage, related to the symmetries of the system. The reason is that the action is a scalar. I.e. the degrees of freedom are combined in the action such that their non-trivial behavior under space-time symmetries exactly cancels, leaving the action invariant. In the presence of internal symmetries, the d.o.f. are combined such that the action also remains invariant under these symmetry transformations. These invariance properties are not straight forward to ascertain on the level of the governing equations, where an individual degree of freedom may change in a non-trivial fashion under space-time or internal symmetry transformations. One pertinent example are electric and magnetic fields, which transform into each other when changing between inertial frames with different relative velocities (Lorentz boosts). In contrast, the classical action of electromagnetism remains invariant under all types of Lorentz transformations.

In this study we show that the action formulation offers additional benefits also when it comes to discretization, which is the main reason why we set out to formulate and solve classical IBVPs directly on the level of the system action.

It is important to note that the action from introductory classical field theory is derived in the context of a boundary value problem, where the field configuration is known at both the initial and final times. Such a BVP setting is non-causal, as the final state must be known before the intermediate field configurations can be computed. In order to exploit the benefits of the action formalism for genuine initial boundary value problems, where only initial and boundary data is known, it has been shown in [14] and put in practice in [15] that one must construct an action with a doubled number of degrees of freedom. (This construction will be discussed in more detail in <ref> and <ref>)

As described in the opening motivation, our aim is to address three central challenges associated with the discrete treatment of IBVPs: 1) loss of space-time symmetries, 2) appropriate choices of meshes and 3) the construction of boundary conditions.

The breaking of translation symmetry is one key artifact in discretizing space-time and representative of the first challenge. In the presence of a finite grid spacing, only translations of that size (or integer multiples thereof) can be implemented, prohibiting arbitrary infinitesimal translations. In turn the concept of a continuous symmetry, the basis for Noether's theorem [16], is absent and the associated Noether charges, i.e. energy, linear and angular momentum, of the system are no longer conserved.

In some systems it may be possible to partially salvage the situation with regard to time translations, by going over to a Hamiltonian formulation, where time at first is kept continuous. But the symplectic solvers deployed in those cases only conserve energy on average (see e.g. [17]). In addition, the Hamiltonian approach relies crucially on the availability of well-defined canonical momenta. Especially in systems with inherent constraints (see the discussion in [18]), these may be difficult to define or may require the choice of a particular gauge, as is the case with the Gauss constraint in classical electrodynamics.

Energy conservation is intimately tied to the stability of the system dynamics (see e.g. [19] and also [20]), hence non-conservation adversely affect the stability of numerical schemes. There is research ongoing to define a Noether theorem for finite grids in the context of so-called calculus of variations on time scales (see e.g. [21]). It however requires one to introduce a specific discretization scheme, which does not possess the same beneficial mimetic[A mimetic discretization exactly mimics certain continuum properties in the discrete setting.] properties as e.g. summation-by-parts (SBP) operators (for reviews see e.g. [22, 23, 24]), and which appears difficult to extend to higher orders. To give Noether's theorem inherent meaning requires realizing a discretization in which infinitesimal space-time symmetry transformations can be accommodated.

The second challenge is related to the appropriate construction of meshes on which the solution of the discretized IBVP is obtained. It is intuitively clear that regions in space-time where the field evolution is more rapid require a higher density of grid points to maintain accuracy of the solution. The systematic construction of these grids is known as adaptive mesh refinement (AMR)  [25, 26, 27], a well established field of research with a rich history. Over the past decades different strategies for the refinement of meshes have been put forward, based on three central concepts. The first is feature detection (see e.g. [28]), which deploys so-called sensors (often for shocks) to determine where violent dynamics require increased resolution. The second follows the evolution of so-called adjoint variables (see e.g. [29, 30]) that can be used to improve the accuracy of the solution itself, as well as to provide an error bound, informing one in which space-time region accuracy of the solution is lost. The third approach is based on aposteriori error estimates (see e.g. [31, 32, 33]) which can be used to direct refinement of meshes to maintain a desired error threshold.

Each of the above strategies thus proposes a different set of criteria to identify space-time regions in which the resolution of the simulation must be modified, given a certain user defined goal for accuracy. It would be highly desirable to uncover a theoretical guiding principle behind mesh refinement, which  produces an inherent AMR criterion directly from the formulation of the IBVP.

The third challenge arises from difficulties in constructing boundary conditions for the simulation domain. In particular in the handling of propagating features, such as localized wave packets or extended wave fronts, unphysical reflections at the boundaries may and most often do contaminate the interior of the domain. Several approaches to so-called non-reflecting boundary conditions have been developed, among them complete radiation boundary conditions (CRBC) [34, 35] and perfectly matched layers (PML) [36].
Partially non-reflecting boundaries may be constructed, such that they allow wave-fronts traveling in a certain direction to pass. While it is already possible to treat systems of waves with different dispersion relations (see e.g. [37]), so far no conclusive solution to the challenge of non-reflecting boundaries has been obtained. In order to support a comprehensive resolution of this challenge, formulations that provide additional flexibility of boundary treatment are called for.

In order to address the first two of the three above mentioned challenges, our study puts forward two central contributions.

The first central contribution of this work is the construction of a novel continuum action for second order IBVPs in $(d+1)$ dimensions. Instead of being formulated directly in space-time, we introduce a $(d+1)$ dimensional abstract parameter space in which all degrees of freedom, including space-time coordinate maps, propagate. These coordinate maps translate parameter coordinates into space-time coordinates and enter the action together with the propagating fields.

The presence of dynamical coordinate maps leads to an interplay with the fields, which determines how fast space-time progresses in terms of the abstract coordinates. The fact that coordinate maps are now dynamical means that they too appear in the boundary terms one encounters, when deriving the governing equations from the action. Their presence thus offers new flexibility in how to implement consistent boundary conditions and thus offers new possibilities for developing non-reflecting boundary conditions.

The second central contribution is a strategy to discretize the novel action in the abstract parameters and not in space-time coordinates. To this end we deploy appropriately regularized multidimensional summation-by-parts finite difference operators. This discretization takes place on the level of the abstract parameters with a naive hypercubic mesh[The discretization performed here constitutes a multidimensional generalization of the discretization strategy developed for the world-line approach in [38].]. Since the values of the coordinate maps remain continuous, the discretized action also retains its invariance under all infinitesimal continuum space-time symmetries. Explicitly preserving the space-time symmetries in this way ensures that Noether's theorem remains fully valid in the discrete setting and that the associated Noether charges remain conserved.

The interplay between fields and coordinate maps in the novel action leads to a self-consistent regulation of how space-time flows. This furnishes the basis for achieving automatic AMR after discretization. We will demonstrate in <ref>, using scalar wave propagation in $(1+1)$ dimensions as example, that a discretization on the level of the abstract parameter space leads to a non-trivial evolution of the coordinate maps, coupled to the field dynamics. This in turn defines a non-trivial mesh of space-time coordinates. It turns out that our procedure leads to coordinate maps that evolve such that space-time is resolved more finely in regions of rapid dynamics, while space-time resolution is coarser where the dynamics of the fields are less rapid.

The rest of the paper is structured as follows. In <ref> we derive the novel reparameterization invariant continuum action for fields in $(d+1)$ dimensions. We start out in <ref> in the context of the conventional boundary value setting of classical field theory, to highlight the main novelty of our action, which is the presence of dynamical coordinate maps. In <ref> we proceed to formulate the novel action with doubled degrees of freedom, necessary to describe a genuine causal IBVP. In preparation for discretizing the action, we discuss in <ref> how initial and boundary data can be explicitly included in the action via Lagrange multipliers and how this affects the Noether charges of the system. <Ref> lays out the details of our discretization strategy for the novel action based on appropriately regularized summation-by-parts finite difference operators. All theoretical development in the preceding sections will be put to work in <ref>, addressing the concrete example of scalar wave propagation in $(1+1)$ dimensions. We provide the explicit form of the continuum action, discretize this action and solve for the classical solution by numerically locating the critical point of the IBVP action, i.e. without solving the governing equations. The convergence properties of the discretization strategy will be assessed.

§ CONTINUUM ACTION WITH DYNAMICAL COORDINATE MAPS

In this section we construct a reparameterization invariant action for classical field theory with dynamical coordinate maps, inspired by the world-line formalism of relativistic point mechanics (see, e.g. [6, 9]).

For the world-line formalism, the theory of general relativity and its guiding principle of reparameterization invariance provides us with an action that naturally contains dynamical coordinate maps.
Based on this starting point, it is possible to explore (as was done in [38]) the consequences of the presence of coordinate maps, which even after discretization allows for preservation of space-time symmetries. While in string theory actions with explicit coordinate maps have been studied for a long time, e.g. in the context of Born-Infeld electromagnetism on D-branes [39], field theory actions are predominantly formulated with reference to space-time coordinates. Our goal in this section is thus to reverse engineer an action for a wide variety of IBVPs in the presence of dynamical coordinate maps, using the world-line formalism as guide. This constitutes a genuinely novel approach to IBVPs not explored in the literature so far.

Conventionally, the variational principle in classical mechanics and field theory is formulated as a boundary value problem (BVP). One assumes that the values of the propagating degrees of freedom are known at initial time and final time and the critical point of the classical action provides us with the classical trajectory or field configuration spanning between these time slices. As discussed in [14] and more recently in [15], such a BVP is not causal, as the final state of the system must be known before the trajectory can be determined. In order to solve a genuine initial boundary value problem, a different variational treatment is called for, which involves a doubling of the degrees of freedom of the system.

In the following subsection we begin the construction of our novel action with dynamical coordinate maps in the context of classical field theory as a boundary value problem. This allows us to present to the reader the novel ingredients in our approach, before implementing a genuine IBVP formulation in the subsequent subsection.

§.§ Boundary Value Problem formulation

The starting point for our construction of a reparametrization invariant action with dynamical coordinate maps is an analogy with the relativistic dynamics of point particles in general relativity.

Let us briefly review the motion of a point particle under the influence of arbitrary forces, encoded in the potential function $V(x)$ (see [15] for a more detailed discussion). The standard non-relativistic action, applicable for motion at velocities much slower than the speed of light reads

\begin{align}
S_{\rm nr}=\int \, dt\, L_{\rm nr}[\vec{x},\dot{\vec{x}}] = \int \, dt\, \Big\{ \frac{1}{2} m \dot \vec{x}^2(t) - V\big(\vec{x}(t)\big) \Big\}, \label{eq:nonrelac}
\end{align}

where $L_{\rm nr}$ refers to the classical Lagrangian, which consists of the difference between kinetic and potential energy.
In contrast, the world-line formalism of general relativity posits that in order to correctly describe the physics also at velocities close to those at the speed of light, one must consider the following action

\begin{align}
S_{\rm wl}&=\int \, d\gamma\,\Big\{  (-mc)\sqrt{ G_{\mu\nu}\frac{d X^\mu}{d\gamma} \frac{dX^\nu}{d\gamma} } - V(X) \Big\}. \label{eq:geodesic}
\end{align}

Summation over repeated indices is implied. The kinetic term of this action is derived following the guiding principle of reparameterization invariance, central to the general theory of relativity[In order to show that the potential term too obeys reparameterization invariance, one has to start from the so-called Polyakov action with an auxiliary field $\xi$, which allows one to couple mass and potential terms in a manifestly reparameterization invariant fashion. The Polyakov action can be shown to be classically equivalent to the world-line action [40].]. The action $ S_{\rm wl}$ describes the motion of a point particle in a space-time characterized by the metric tensor $G$[The metric tensor $G$ describes how length and time intervals, as well as angles are defined in a given space-time. In flat $(3+1)$ dimensional space-time $G\in\mathbb{R}^{4 \times 4}$ is often denoted as $\eta_{\mu\nu}\equiv G_{\mu\nu}$, has negative determinant (we choose here the signature convention $(+,-,-,-)$ [10]) and is diagonal with the entries $G_{\mu\nu}={\rm diag}[G_{00},G_{11},G_{22},G_{33}]={\rm diag}[c^2,-1,-1,-1]$. The inverse of $G_{\mu\nu}$ is denoted by $G^{\mu\nu}$.], but instead of time $t$, uses an abstract world-line parameter $\gamma$ to parameterize the one-dimensional trajectory $X^\mu(\gamma)=(t(\gamma),\vec{x}(\gamma))^\mu=(t(\gamma),x_1(\gamma),\ldots,x_d(\gamma))^\mu$ that the particle traces out in $(d+1)$ dimensional space-time. Here the superscript $\mu=0,\ldots,d$ refers to the $\mu$-th components of the vector $X$ and is used in the following to implement the Einstein summation convention. As discussed in [15], for a diagonal metric and a potential, which is small compared to the rest energy of the point particle $V(\vec{x})/2mc^2\ll1$, one may rewrite the world-line action <ref> as a genuine geometric, i.e. a geodesic action

\begin{align}
S_{\rm wl}\approx S_{\rm geo} = \int \, d\gamma\,(-mc)\Big\{  \sqrt{ \Big(G_{00}+\frac{V(\vec{x})}{2mc^2}\Big)\frac{d X^0}{d\gamma} \frac{dX^0}{d\gamma}+ G_{ii}\frac{d X^i}{d\gamma} \frac{dX^i}{d\gamma} } \Big\}.\label{eq:geodesicincl}
\end{align}

When taking the non-relativistic limit in flat space-time $G=\eta$ (see <ref>), which builds on the separation of two scales $V(\vec{x})/2mc^2\ll 1$ and $\frac{d|\vec{x}|/d\gamma}{dt/d\gamma}/c = v/c \ll 1$, one is able to recover a non-relativistic action from <ref>. This non-relativistic action features one decisive difference with respect to the usual <ref>

\begin{align}
S_{\rm nr} \overset{\rm from \; \cref{eq:geodesicincl}}{=} \int \, dt\, \Big\{-mc^2 +  \frac{1}{2} m \dot{\vec{x}}^2(t) - V\big(\vec{x}(t)\big) \Big\},\label{eq:nonrelacderiv}
\end{align}

which is the presence of the constant rest mass term $-mc^2$. In addition one also learns that the factor $1/2$ in front of $m\dot{\vec{x}}^2$ is directly related to the expansion of the square root in the world-line action.

Note that a constant term in the action does not affect the classical trajectory. Indeed, a constant term vanishes under variation of the action, and does not contribute to the stationary condition $\delta S_{\rm nr}=0$ defining the trajectory. This explains why one cannot guess the “correct value" of the constant term without an additional guiding principle. In the case of general relativity, the guiding principle is reparametrization invariance, which leads to the construction of the world-line action, and the non-relativistic limit $v\ll c$ is a systematic approximation of that world-line action. In particular, this expansion in $v/c$ leads to a series of correction terms to the non-relativistic action <ref>. Since the effect of these corrections only becomes relevant in the genuine relativistic regime, it was historically difficult in experiment to spot that <ref> is missing relevant physics.

Keeping the above line of arguments in mind, let us turn to genuine field theory. The fully relativistic field theory action found in standard textbooks (see e.g. chapter 4 of [9]) includes both kinetic terms, referring to the space-time derivatives of the field, as well as (self-)interaction terms in the form of an in general non-linear potential function $V(\phi)$. This action too is reparameterization invariant and is conventionally written as

\begin{align}
S=\int d^{(d+1)}X \, \sqrt{-{\rm det}[G]} \frac{1}{2}\Big( G^{\mu\nu} \partial_\mu\phi(X)\partial_\nu\phi(X) - V(\phi) \Big). \label{eq:stdFTaction}
\end{align}

For notational convenience, we took the liberty to define the self-interaction function $V(\phi)$ such that an overall factor $1/2$ is taken to the front (i.e. a mass term quadratic in the fields would now read $V_{\rm mass}=m\phi^2$). The above action describes second order systems, which feature a kinetic term quadratic in the velocities $\partial_\mu\phi$, as is the case for e.g. scalar wave propagation [4] (for which $V(\phi)=0$) or self interacting Klein-Gordon fields with $V(\phi)\neq0$.

To connect the action of <ref>, written in an abstract relativistic form, to the familiar physics of, e.g., scalar wave propagation, let us consider a flat space-time, where ${\rm det}[G]=-c^2$, and let us set $V(\phi)=0$. We then recover the standard action for scalar wave propagation

\begin{align}
S^{\rm wp}=\int d^{(d+1)}X\, c\,\underbracket{\Big( \frac{1}{c^2}\,\frac{1}{2}  \partial_0\phi(X)\partial_0\phi(X) - \frac{1}{2}\partial_i\phi(X)\partial_i\phi(X)  \Big)}_{{\cal L}^{\rm wp}}. \label{eq:stdWPaction}
\end{align}

The factor $c$ up front in the integrand ensures that the units of the integration measure match those of the derivatives in the parentheses. The variational principle tells us that the classical solution is obtained at the critical point of this action. After computing the stationary condition

\begin{align}
&0=\delta S^{\rm wp}= \int d^{(d+1)}Xc \Big\{ \frac{1}{c^2} (\partial_0\phi)\delta(\partial_0\phi) - (\partial_i\phi)\delta(\partial_i\phi)\Big\}\\
\nonumber &\overset{\rm IBP}{=}\int d^{(d+1)}Xc \Big\{ -\Big( \frac{1}{c^2}  \partial^2_0\phi(X) - \partial^2_i\phi(X)  \Big)\delta\phi \Big\} + \frac{1}{c^2}(\partial_0\phi \,\delta\phi)|_{t^{\rm i}}^{t^{\rm f}} + \sum_j  (\partial_j\phi\, \delta \phi)|_{x_j^{\rm i}}^{x_j^{\rm f}}
\end{align}

one finds that only if the spatial boundary terms vanish (the temporal ones in the BVP setting do, since $\delta \phi$ vanishes), one is led to the following governing (Euler-Lagrange) equation

\begin{align}
\frac{\partial {\cal L}^{\rm wp}}{\partial  \phi(X)} - \partial_\mu \Big[ \frac{\partial {\cal L}^{\rm wp}}{\partial (\partial_\mu \phi(X))}\Big]= -\Big( \frac{1}{c^2}  \partial^2_0\phi(X) - \partial^2_i\phi(X)  \Big)=0, \label{eq:stdWPgoveq}
\end{align}

i.e. the well-known scalar wave equation with wave propagation speed $c$.

The standard action in <ref> remains invariant under arbitrary differentiable reparameterizations of the space-time $X\to X^\prime$ due to the presence of the term $\sqrt{-{\rm det}[G]}$. This term absorbs the Jacobian $J=\partial X/\partial X^\prime$ of the coordinate transformation to yield $\sqrt{-{\rm det}[J^TGJ]}=\sqrt{-{\rm det}[G^\prime]}$, which defines the new induced metric $G^\prime$ in the new coordinate system[The metric $G$ is associated with the inner product of four-dimensional space-time vectors, which in turn defines distances between events taking place in the space-time described by the coordinates $X$. The induced metric $G^\prime$ implements the corresponding inner product in the coordinate system defined by $X^\prime$. For more details see, e.g., [40].]. The action is also invariant under the full Poincare group of space and time translations, spatial rotations and relativistic rotations involving space and time, called boosts.

After inspecting the properties of the standard action, let us take the first step towards reverse engineering a field theory action with dynamical coordinate maps. Taking insight from the relation between the non-relativistic point-particle action <ref> and the geodesic action <ref>, we add to <ref> a constant term $-T$ akin to the rest mass term $-mc^2$ (or tension $T$ of a D-brane, see discussion on page par:tension) . Since the extremum of $S_{\rm ft}$ is unchanged by the constant, the classical physics too remains unchanged

\begin{align}
S=\int d^{(d+1)}X \, \sqrt{-{\rm det}[G]} \Big\{ -T + \frac{1}{2}\Big(G^{\mu\nu} \partial_\mu\phi(X)\partial_\nu\phi(X) + V(\phi) \Big)\Big\}.\label{eq:stdFTactionpconst}
\end{align}

Note that in order for the constant $T$ to be accommodated in this action, it must carry the units of an energy density.

For the next step, we interpret the presence of the constant term and the factor $1/2$ as an indication that we are actually dealing with the low energy limit of a more general action (akin to the world-line action) which contains a square root term. I.e. we assume that in <ref> small terms at higher order in inverse powers of some dimensionless ratio $\kappa\sim{\rm energy\, density}/T$ have been discarded. As part of our reverse engineering strategy, we undo this expansion by introducing a square root as

\begin{align}
\nonumber&{\cal S}_{\rm BVP}\equiv\\
&\int d^{(d+1)}X \, \sqrt{-{\rm det}[G]} \big(-T\big) \Big\{ 1 - \frac{1}{2T}\Big(G^{\mu\nu} \partial_\mu\phi(X)\partial_\nu\phi(X) + V(\phi) \Big)\Big\} + {\cal O}(\kappa^2),\\
&=\int d^{(d+1)}X \, \sqrt{-{\rm det}[G]} \big(-T\big) \sqrt{ 1 - \frac{1}{T}\Big(G^{\mu\nu} \partial_\mu\phi(X)\partial_\nu\phi(X) + V(\phi) \Big)}.  \label{eq:stdFTactionbacksqrt}
\end{align}

In order for the previous step to be justified the ratio $\kappa$ must indeed be small; conversely, $T$ must correspond to a large energy density compared to the energy density represented by the field $\phi$. This situation mirrors the world-line formalism, where the rest-mass of a point particle is much larger than the kinetic and potential energy encountered on everyday scales.

Now we take the next vital step, which consists of elevating the space-time coordinates $X^\mu$ to dynamical degrees of freedom. This too mirrors the world-line formalism, where one elevates $t\to t(\gamma)$ and it provides the mechanism that will subsequently enable the preservation of space-time symmetries after discretization.

To this end we change coordinates in the integral from space-time $X^\mu=(t, \vec{x})^\mu=(t,x_1,\ldots,x_d)^\mu$ (in <ref>) to a set of arbitrary parameters $\Sigma^a=(\tau,\vec{\sigma})^a=(\tau,\sigma_1,\ldots,\sigma_d)^a$ with $a=0,\ldots,d$ (in <ref>). The two coordinate systems are related by an invertible and differentiable mapping $X^\mu(\Sigma)$, which is explicitly encoded via the Jacobian $J^\mu_a=\partial X^\mu/\partial \Sigma^a$ that appears in the measure as we carry out the change of variables

\begin{align}
{\cal S}_{\rm BVP}=\int d^{(d+1)}\Sigma \, |{\rm det}[J]|&\sqrt{-{\rm det}[G]} \big(-T\big)\label{eq:stdFTactionbacksqrt2}\\
\nonumber&\times \sqrt{ 1 - \frac{1}{T}\Big(G^{\mu\nu} \partial_\mu\phi(X(\Sigma))\partial_\nu\phi(X(\Sigma)) + V(\phi) \Big)}.
\end{align}

We will take the range of the parameters $\Sigma^a$ to be finite and to span a hypercube interval from their lowest $(\tau^{\rm i},\sigma_1^{\rm i},\ldots,\sigma_d^{\rm i})$ to highest values $(\tau^{\rm f},\sigma_1^{\rm f},\ldots,\sigma_d^{\rm f})$. I.e. we will always refer to initial and final values of the parameters with superscript letters ${\rm i}$ and ${\rm f}$ respectively.

Let us inspect the induced metric $g_{ab}$ on the parameter manifold, defined as the quantity $g_{ab}\equiv G_{\mu\nu} J^\mu_a J^\nu_b$. The induced metric naturally arises from combining the determinant of the Jacobian with the metric dependent contribution of the measure

\begin{align}
&\sqrt{-{\rm det}[J]{\rm det}[G]{\rm det}[J]}\\
\nonumber &\qquad\qquad=\sqrt{-{\rm det}[J^T]{\rm det}[G]{\rm det}[J]}=\sqrt{-{\rm det}[J^T G J]}= \sqrt{-{\rm det}[g]},
\end{align}

so that the action reads

\begin{align}
{\cal S}_{\rm BVP}=\int d^{(d+1)}\Sigma \, \sqrt{-{\rm det}[g]}& \big(-T\big)\label{eq:stdFTactionbacksqrtmod1}\\
\nonumber&\times \sqrt{ 1 - \frac{1}{T}\Big(G^{\mu\nu} \partial_\mu\phi(X(\Sigma))\partial_\nu\phi(X(\Sigma)) + V(\phi) \Big)}.
\end{align}

If we also absorb the determinant ${\rm det}[g]$ into the main square root expression in the action, we obtain

\begin{align}
{\cal S}_{\rm BVP}&=\int d^{(d+1)}\Sigma \, \big(-T\big) \sqrt{ -{\rm det}[g] + \frac{1}{T}\Big(G^{\mu\nu} \partial_\mu\phi(\Sigma)\partial_\nu\phi(\Sigma) + V(\phi) \Big){\rm det}[g]}.  \label{eq:stdFTactionbacksqrtmod2}
\end{align}

What we have achieved so far is to formulate an action for fields (<ref>) in an arbitrary parameter space, described by the coordinates $\Sigma^a$. The maps $X(\Sigma)$ that relate coordinates in parameter space to space-time coordinates enter in the ${\rm det}[g]$ terms in the square root. Note that the field dependent term on the right in <ref> however still contains reference to the metric $G$ and derivatives in terms of the old coordinates $\partial_\mu=\partial/\partial X^\mu$. In order to complete the transition from $X$ to $\Sigma$ coordinates, we have to rewrite the RHS of <ref>  solely in terms of derivatives in the new coordinates $\partial_a=\partial/\partial \Sigma^a$.

To this end let us recognize that the contraction of the partial derivatives $G^{\mu\nu}\partial_\mu\partial_\nu$ is invariant under reparameterization. It thus follows that we can simply replace it by a contraction in terms of the new derivatives  $\partial_\mu\phi \partial_\nu\phi G^{\mu\nu}=\partial_a\phi \partial_b\phi g^{ab}$. We can confirm this replacement by explicit calculation, writing in component notation $(\partial^a\phi J_a^\mu) (\partial^b\phi J_b^\nu) G_{\mu\nu}$ with $J^a_\mu J^b_\nu G^{\mu\nu}=g^{ab}$ by definition of $g$ and the fact that $J^a_\mu=(J^{-1})^\mu_a$. Let us also move the potential term to the left under the square root, then

\begin{align}
{\cal S}_{\rm BVP}&=\int d^{(d+1)}\Sigma \, \big(-T\big) \sqrt{ \Big(\frac{1}{T}V(\phi)-1\Big){\rm det}[g] + \frac{1}{T} \partial_a\phi(\Sigma)\partial_b\phi(\Sigma) g^{ab}  {\rm det}[g]}.  \label{eq:stdFTactionbacksqrtmod3}
\end{align}

Note that both the inverse (denoted by $g^{-1}=g^{ab}$) and the determinant ${\rm det}[g]$ of the induced metric $g_{ab}$ appear in the $\phi$ dependent term on the right under the square root. As long as $g$ is invertible, we can further simplify the action by applying a well-known matrix identity, involving the adjugate of the metric, ${\rm adj}[g]=g^{-1}{\rm det}[g]$. This leads us to the final expression for our novel reparameterization invariant action with dynamical coordinate maps

\begin{align}
\hspace{-0.8cm}\tcbhighmath[boxrule=2pt,arc=1pt,colback=gray!10!white,colframe=gray,
drop fuzzy shadow=gray]{{\cal S}_{\rm BVP}=\int d^{(d+1)}\Sigma \, \big(-T\big) \sqrt{ \Big(\frac{1}{T}V(\phi)-1\Big){\rm det}[g] + \frac{1}{T} \partial_a\phi(\Sigma)\partial_b\phi(\Sigma) {\rm adj}[g]_{ab}}.} \label{eq:novelaction}
\end{align}

Summation of the indices $a$ and $b$ is implied. Since no additional factors of $g$ enter this summation, we write both indices as subscripts. We emphasize that the induced metric $g$ makes reference to $X^\mu$ and $\partial_aX^\mu$. The corresponding Lagrangian ${\cal L}_{\rm BVP}$ is defined as ${\cal S}_{\rm BVP}\equiv\int d^{(d+1)}\Sigma\,{\cal L}_{\rm BVP}[X,\partial_a X,\phi,\partial_a \phi]$.

This novel action <ref> constitutes a generalisation of the world-line action for a point particle to field dynamics. Thus we may expect some similarities but also differences between <ref> and <ref>.

Note first that our new action is manifestly invariant under Poincaré transformations. These transformations only affect the coordinate maps $X^\mu$, which appear in the form of the determinant or adjugate of the matrix $g$. As the coordinate maps enter the induced metric $g$ by contraction of their derivatives via the metric tensor $G$, the entries of $g$ by themselves are invariant under Lorentz transformations. And since it is only the derivatives of $X^\mu$ that appear, any constant shift in the coordinate maps is annihilated. In turn the entries of $g$ and hence ${\cal S}$ are invariant under the full Poincaré group of transformations.

Note also the overall multiplicative factor of $-T$ in <ref>, the analog to $-mc$ for the world-line formalism. This quantity in both cases denotes the scale at which the full dynamical nature of the coordinates becomes relevant. For a point particle, at kinetic energies close to its rest energy, we must take into account that its motion through time can no longer be considered independent from its motion through space, contrary to the non-relativistic limit. Similarly here, once the energy density carried by the field becomes of the order of $T$, the dynamics of the coordinate mappings $X^\mu$ encoded in $g$ are inescapably intertwined with those of the field $\phi$.

In the geodesic action <ref>, the effects of a potential were approximately absorbed into a term related to the dynamics of the coordinate $t$, promoted to a dynamical degree of freedom. In the case of field theory in <ref>, the potential function $V(\phi)$ is naturally associated with the ${\rm det}[g]$ term under the square root, which governs the newly introduced dynamical coordinate maps.

Let us consider the physics of the coordinate maps in the absence of a field, by taking the $T\to\infty$ limit. Interestingly, in that case, we recover an action given solely by the square root of the determinant of the induced metric. When $d=1$, this expression plays a central role in string theory, where it is known as the Nambu-Goto action (for details see [40]) and describes a spatially extended string that propagates, tracing out a so-called world-sheet. In an otherwise flat space-time, the string world-sheet, as the name suggest, constitutes a flat sheet with constant slope. For $d>1$, this expression is known as a D-brane action (see, e.g., [41]).

In our case the action includes reference to all $(d+1)$ space-time coordinates $X$ and we consider a $(d+1)$-dimensional world-volume that is traced out, as the underlying parameters $\Sigma$ vary. The constant $T$ within string theory has a clear physics interpretation: it denotes the tension $T/c$ of the string, and constitutes a fundamental property of nature. Similarly we may interpret $T$ here as a reference energy density that tells us how much energy is required in order to produce a stretching or compression of the coordinate mappings.

The D-brane action is invariant under coordinate reparameterizations and by construction, so is our novel action, including the field degrees of freedom.

The variational principle of classical field theory tells us that the classical field configuration is found as the critical point of the system action. Furthermore this critical point can either be found by direct variation of the action itself or equivalently by solving a set of partial differential equations of motion, the so-called Euler-Lagrange equations.

To obtain the Euler-Lagrange equations of motion, let us proceed by carrying out the variation of <ref> in all the dependent variables $X^\mu$, $\partial_a X^\mu$, $\phi$ and $\partial_a \phi$.

\begin{align}
\nonumber&\delta {\cal S}_{\rm BVP}[X, \partial_a X, \phi, \partial_a \phi]\\
\nonumber &=\int d^{(d+1)}\Sigma \Big\{
\frac{\partial {\cal L}_{\rm BVP}}{\partial X^\mu}\delta X^\mu + \frac{\partial {\cal L}_{\rm BVP}}{\partial (\partial_a X^\mu)}\delta (\partial_a X^\mu) + \frac{\partial {\cal L}_{\rm BVP}}{\partial \phi}\delta \phi + \frac{\partial {\cal L}_{\rm BVP}}{\partial (\partial_a \phi)}\delta (\partial_a \phi) \Big\}\\
&=\int d^{(d+1)}\Sigma \Big\{
\Big( \frac{\partial {\cal L}_{\rm BVP}}{\partial X^\mu} -\partial_a \frac{\partial {\cal L}_{\rm BVP}}{\partial (\partial_a X^\mu)}\Big)\delta X^\mu + \Big(  \frac{\partial {\cal L}_{\rm BVP}}{\partial \phi} - \partial_a \frac{\partial {\cal L}_{\rm BVP}}{\partial (\partial_a \phi)} \Big)\delta \phi \Big\} \label{eq:varSeom}\\
&\qquad + \sum_{n=0}^{d} \int \prod_{\tiny\begin{array}{c}a=0\\a\neq n\end{array}}^{d}d\Sigma_a \Big[ \frac{\partial {\cal L}_{\rm BVP}}{\partial (\partial_n X^\mu)} \delta X^\mu \Big]_{\Sigma_n^{\rm i}}^{\Sigma_n^{\rm f}}\label{eq:varSbnd1}\\
&\qquad + \sum_{n=0}^{d} \int \prod_{\tiny\begin{array}{c}a=0\\a\neq n\end{array}}^{d}d\Sigma_a \Big[ \frac{\partial {\cal L}_{\rm BVP}}{\partial (\partial_n \phi)} \delta \phi \Big]_{\Sigma_n^{\rm i}}^{\Sigma_n^{\rm f}}
\label{eq:varSbnd2}\end{align}

The equivalence between the critical point, defined by $\delta {\cal S}_{\rm BVP}[X, \partial_a X, \phi, \partial_a \phi]=0$ and the equations of motion in the parentheses in (<ref>) requires the boundary terms in  (<ref>) and (<ref>) to vanish. In the BVP setting we consider here, the case $n=0$ makes reference to boundaries in the temporal direction $\Sigma_0=\tau$. By construction, we keep the values of the d.o.f. fixed at initial and final temporal parameter. I.e. the variations $\delta X^\mu$ and $\delta\phi$ have to vanish on the temporal boundaries at initial $\tau^{\rm i}=\Sigma_0^{\rm f}$ and final $\tau^{\rm f}=\Sigma_0^{\rm f}$. For the spatial boundaries related to the lowest and highest values of $\sigma_1,\ldots,\sigma_d$ there exist different possibilities to make the corresponding terms in (<ref>) and (<ref>) vanish. One possibility is to impose Dirichlet boundary conditions, which entails that the variations themselves vanish. On the other hand, another possibility is to impose von-Neumann or mixed boundary conditions, which involve specification of the derivatives of the degrees of freedom. We will discuss the treatment of boundary conditions in more detail in the explicit example of wave propagation in 1+1 dimensions in <ref>.

In preparation for the next steps, recall that the location of the critical point of a functional does not change if we apply a monotonic function to its integrand. Also, for a reparameterization invariant action there exist an infinite number of ways to express the classical dynamics of the system it describes. In the world-line formalism, which describes the motion of a single point particle along its world-line trajectory, reparameterization invariance refers to the fact that the same world-line trajectory in $(d+1)$ dimensions can be traversed with different speeds relative to the world-line parameter $\gamma$. Such freedom leads to difficulties for a numerical minimizer in locating the classical field configuration. At the same time, the invariance of the physics under different parameterizations means that we are free to choose a parameterization at will.

These difficulties were avoided in the world-line formalism by squaring the integrand of the action. Motivated by this choice, we consider an action ${\cal E}$, where the integrand has been squared and the pre-factors have been modified for notational convenience

\begin{align}
\nonumber{\cal E}_{\rm BVP}&= \int d^{(d+1)}\Sigma\; E_{\rm BVP}[X,\partial_a X,\phi, \partial_a\phi]\\
&=\int d^{(d+1)}\Sigma \; \frac{1}{2} \Big\{ \Big(\frac{1}{T}V(\phi)-1\Big){\rm det}[g] + \frac{1}{T} \partial_a\phi(\Sigma)\partial_b\phi(\Sigma) {\rm adj}[g]_{ab}\Big\}. \label{eq:novelactionE}
\end{align}

Note that this action, as intended, is no longer invariant under general reparameterizations of $\Sigma$. However the critical point of ${\cal E}_{\rm BVP}$ and ${\cal S}_{\rm BVP}$ are equivalent, as the integrand in the action differs by only a monotonic transformation.

The equations of motion associated with the functional are obtained from its variation as

\begin{align}
\nonumber &\delta {\cal E}_{\rm BVP}[X, \partial_a X, \phi, \partial_a \phi]\\
&=\int d^{(d+1)}\Sigma \Big\{
\frac{\partial E_{\rm BVP}}{\partial X^\mu}\delta X^\mu + \frac{\partial E_{\rm BVP}}{\partial (\partial_a X^\mu)}\delta (\partial_a X^\mu) + \frac{\partial E_{\rm BVP}}{\partial \phi}\delta \phi + \frac{\partial E_{\rm BVP}}{\partial (\partial_a \phi)}\delta (\partial_a \phi) \Big\}\label{eq:varE}\\
&=\int d^{(d+1)}\Sigma \Big\{
\Big( \frac{\partial E_{\rm BVP}}{\partial X^\mu} -\partial_a \frac{\partial E_{\rm BVP}}{\partial (\partial_a X^\mu)}\Big)\delta X^\mu + \Big(  \frac{\partial E_{\rm BVP}}{\partial \phi} - \partial_a \frac{\partial E_{\rm BVP}}{\partial (\partial_a \phi)} \Big)\delta \phi \Big\} \label{eq:varEeom}\\
&\qquad + \sum_{n=0}^{d} \int \prod_{\tiny\begin{array}{c}a=0\\a\neq n\end{array}}^{d}d\Sigma_a \Big[ \frac{\partial E_{\rm BVP}}{\partial (\partial_n X^\mu)} \delta X^\mu \Big]_{\Sigma_n^{\rm i}}^{\Sigma_n^{\rm f}}\label{eq:varEbnd1}\\
&\qquad + \sum_{n=0}^{d} \int \prod_{\tiny\begin{array}{c}a=0\\a\neq n\end{array}}^{d}d\Sigma_a \Big[ \frac{\partial E_{\rm BVP}}{\partial (\partial_n \phi)} \delta \phi \Big]_{\Sigma_n^{\rm i}}^{\Sigma_n^{\rm f}}
\label{eq:varEbnd2}\end{align}

under the condition that the boundary terms (<ref>) and (<ref>) vanish.

For a manifestly Poincaré invariant action, such as <ref>, Noether's theorem guarantees that each of its global continuous symmetries leads to one conserved current, which in turn defines a conserved charge. We may use the fact that $\delta_X {\cal E}\equiv{\cal E}[X+\delta X]-{\cal E}[X]=0$ and since the transformation only affects changes in the coordinate maps one ends up with the first two terms in the curly brackets in (<ref>). Carrying out integration by parts but retaining the total derivative we obtain

\begin{align}
0=\delta_X {\cal E}_{\rm BVP}&=\int d^{(d+1)}\Sigma \Big\{
\underbracket{\Big( \frac{\partial E_{\rm BVP}} {\partial X^\mu} -\partial_a \frac{\partial E_{\rm BVP}}{\partial (\partial_a X^\mu)}\Big)}_{\rm equations\,of\,motion\, for \, X}\delta X^\mu  +  \partial_a\Big[ \underbracket{\frac{\partial E_{\rm BVP}}{\partial (\partial_a X^\mu)} \delta X^\mu}_{J^a} \Big]\Big\}\label{eq:derivN1}\end{align}

Noether's theorem only holds if the equations of motion are fulfilled, i.e. if the term in the parenthesis in (<ref>) vanishes. We see that the conserved Noether current $J^a$ then is given by the expression

\begin{align}
\partial_aJ^a= \partial_0J^0 + \sum_{n=1}^{d} \partial_n J^n =   \partial_a\Big[ \frac{\partial E_{\rm BVP}}{\partial (\partial_a X^\mu)} \delta X^\mu \Big] =0 \label{eq:NoetherCur}
\end{align}

The corresponding conserved charge is also defined as in regular field theory,

\begin{align}
Q= \int \prod_{a>0}^d d\Sigma^a J^0 = \int d\sigma_1 \ldots d\sigma_d J^0, \label{eq:defQ}
\end{align}

as the spatial volume average over the temporal component of the Noether current $J^0$. Note that reference to temporal and spatial here is understood to refer to the parameters $\tau$ and $\vec{\sigma}$ respectively. In order for $Q$ to be constant in time, its temporal derivative must vanish. To evaluate $\partial_\tau Q$ we first exploit the definition <ref> and then apply the fundamental theorem of calculus

\begin{align}
\partial_\tau Q&\overset{def.\,\cref{eq:defQ}}{=} \int   d\sigma_1 \ldots d\sigma_d \big(\partial_\tau J^0\big) \overset{\cref{eq:NoetherCur}}{=} - \int   d\sigma_1 \ldots d\sigma_d \big( \sum_{n=1}^{d} \partial_{\sigma_n}J^n\big) \\
&=-\sum_{n=1}^{d} \int \prod_{\tiny\begin{array}{c}a=0\\a\neq n\end{array}}^{d}\left.d\sigma_a J^n \right|_{\sigma_n^{\rm i}}^{\sigma_n^{\rm f}}.
\end{align}

One has to inspect how the spatial components $J^n$ behave on the spatial boundaries. Only if they vanish there, does <ref> imply $\partial Q/\partial t=0$ (c.f. (<ref>) and (<ref>)). We will give an explicit expression for the relevant Noether charges when exploring 1+1 dimensional wave motion in <ref>.

§.§ Initial value formulation

The action functional of <ref> does not yet allow us to formulate a causal initial boundary value problem. This becomes apparent, when we consider the derivation of the equations of motion from it in <ref>. The boundary terms only vanish if the variation on the final $\tau$ slice are set to zero, i.e. the values of the d.o.f. must be known at final time.

The way to accommodate the natural absence of information about the final state of the system is to introduce a cleverly arranged doubling of degrees of freedom. This approach to causal initial value problems is known in the quantum physics community as the Schwinger-Keldysh construction (see e.g. [42]) and has been independently derived in the classical context by Galley in [14]. After deploying this construction to conventional point particle motion in [15] and to the world-line formalism in [38], we will now exploit it to formulate a novel action with dynamical coordinate maps for initial boundary value problems involving fields.

The doubling of degrees of freedom can be understood (see sketches and caption in <ref>) to furnish a double shooting method. One considers one copy of the degrees of freedom to live on the so-called forward temporal branch, a second copy on the backward temporal branch. The forward branch is initialized by the physical initial conditions while the backward branch receives information from the final state of the forward branch. As shown and discussed in detail in [14], the solution of the IBVP corresponds to the trajectory, where the forward and backward path agree and the backward path reproduces the correct values and derivatives at initial time.

(left) The conventional BVP formulation of classical field theory, where both initial data $\phi(\tau^{\rm i},\vec{\sigma})$ and final data $\phi(\tau^{\rm f},\vec{\sigma})$ (in green) are provided a priori. While different field configurations exist that connect these two datasets (gray sheets), there exists a unique field configuration $\phi_{\rm cl}$ (in red) that constitutes the critical point of the action. It is this configuration that represents the configuration realized in nature. (right) The doubled d.o.f. construction of [14], necessary to formulate causal IBVPs on the action level. Provided initial data on values $\phi(\tau^{\rm i},\vec{\sigma})$ and derivatives $\dot \phi(\tau^{\rm i},\vec{\sigma})$ (in green), one constructs a double shooting method. One copy of the field evolves forward (arrows to the right) and one copy evolves backward from the final state of the forward branch (arrows to the left). While different field configuration pairs exists that accommodate the initial data on the forward branch (e.g. gray sheets), only a single one (red and blue sheets) fulfills the requirement that it constitutes the critical point of the action with doubled degrees of freedom and that the forward and backward branch solutions agree (for more details see e.g. [15]).

In the following we will thus refer to coordinate mappings and fields on the two branches as $X_1,X_2$ and $\phi_1,\phi_2$, each as functions of the parameters $\Sigma$ (which are not doubled).

Remember that it is the boundary terms in the final $\tau$ slice which prevent a causal formulation in <ref>. The approach of [14] resolves this issue by assigning the forward branch and the backward branch the same Lagrangian with opposite sign in the action

\begin{align}
{\cal E}_{\rm IBVP}&= \int d^{(d+1)}\Sigma\; E_{\rm IBVP}[X_1,X_2,\partial_a X_1,\partial_a X_2,\phi_1, \phi_2, \partial_a\phi_1,\partial_a\phi_2] ,\label{eq:novelactionEIVP}
%\nonumber &=\int d^{(d+1)}\Sigma\; \Big\{ E_{\rm BVP}[X_1,\partial_a X_1,\phi_1, \partial_a\phi_1] -E_{\rm BVP}[X_2,\partial_a X_2,\phi_2, \partial_a\phi_2] \Big\},
\end{align}

\begin{align}
\nonumber E_{\rm IBVP}=E_{\rm BVP}[X_1,\partial_a X_1,\phi_1, \partial_a\phi_1] -E_{\rm BVP}[X_2,\partial_a X_2,\phi_2, \partial_a\phi_2],
\end{align}

and by enforcing the following connecting conditions on the values and derivatives of the degrees of freedom on the final time slice

\begin{align}
\nonumber &X_1^\mu(\tau=\tau^{\rm f},\vec{\sigma})=X_2^\mu(\tau=\tau^{\rm f},\vec{\sigma}), &&\phi_1(\tau=\tau^{\rm f},\vec{\sigma})=\phi_2(\tau=\tau^{\rm f},\vec{\sigma})\\
&\partial_0 X_1^\mu|_{\tau=\tau^{\rm f}}= \partial_0 X_2^\mu|_{\tau=\tau^{\rm f}}, &&  \partial_0 \phi_1|_{\tau=\tau^{\rm f}}= \partial_0 \phi_2|_{\tau=\tau^{\rm f}}.\label{eq:SKconnecting}
\end{align}

As derived in detail in <ref>, these conditions establish a causal variational principle, connecting the equations of motion, with the critical point of the functional ${\cal E}_{\rm IBVP}$ (given that  the spatial boundary terms vanish). The variational principle requires us to provide physical initial conditions about the values and derivatives of the degrees of freedom at initial time as well as connecting conditions between the forward and backward branches, which in turn avoid any non-causal temporal boundary terms.

In setting up the causal variational approach we introduced a second copy of the degrees of freedom. This proliferation of dynamical variables needs to be undone, for which, as argued in [14], one must take the so-called physical limit. This limiting procedure consists of requiring that $X^\mu_1-X^\mu_2=X^\mu_-=0$ and $\phi_1-\phi_2=\phi_-=0$. In turn the dynamic variables $X^\mu_+$, $\phi_+$ converge towards the classical solution. In the equations of motion

\begin{align}
\Big( \frac{\partial E_{\rm IBVP}}{\partial X_\pm^\mu} -\partial_a \frac{\partial E_{\rm IBVP}}{\partial (\partial_a X_\pm^\mu)}\Big)=0, \quad
\Big(  \frac{\partial E_{\rm IBVP}}{\partial \phi_\pm} - \partial_a \frac{\partial E_{\rm IBVP}}{\partial (\partial_a \phi_\pm)} \Big)=0,
\end{align}

the physical limit finds another interpretation. As $E_{\rm IBVP}= E_{\rm BVP}[X_1,\phi_1]- E_{\rm BVP}[X_2,\phi_2]$ arises from the difference of the BVP functionals it contains at least a linear dependence on the relative variables. In the physical limit, thus only the Euler-Lagrange equation where the derivative w.r.t. the relative coordinates is taken, has a chance to survive. Terms in $E_{\rm IBVP}$ that feature at most a linear dependence on $X^\mu_-$ and $\phi_-$ reduce to finite expressions containing $X^\mu_+$ and $\phi_+$. A more rigorous derivation of the physical limit can be made in the Schwinger-Keldysh formalism of quantum field theory, where the physical limit of [14] amounts to taking the classical limit of Feynman's path integral, as discussed in detail in [43].

Let us take a brief look at the type of initial data we must provide to our initial boundary value problem. Since all dynamical quantities are dependent on the parameters $\Sigma^a$ and the kinetic terms in the action <ref> contain derivatives quadratic in $\tau=\Sigma_0$, both the values $X^\mu(\tau^{\rm i},\vec{\sigma}),\phi(\tau^{\rm i},\vec{\sigma})$ and the first derivative in $\tau$, i.e. $\partial_0 X^\mu(\tau^{\rm i},\vec{\sigma}),\partial_0\phi(\tau^{\rm i},\vec{\sigma})$ must be provided. From physical experiments we may obtain access to the values of the field at initial time and we are aware of the space-time coordinates at which this initial field has been measured. We can also determine the temporal derivative $\partial_{X^0}\phi$. On the other hand, the speed by which the coordinate maps and field evolve along $\Sigma_0$ are not fixed a priori. Thus there remains a freedom, as it is only the ratio $\frac{\partial\phi}{\partial X^0} = \frac{\partial \phi}{\partial \Sigma_0}/\frac{\partial X^0}{\partial \Sigma_0}$ that needs to be fixed. Since we consider time $X^0$ as a dynamic variable, its final value at $\tau^{\rm f}$ will depend on the speed assigned to it initially and it is an emergent property of the dynamics between the $X^\mu$'s and $\phi$.

§.§ Lagrange multiplier formulation of initial, boundary and connecting conditions

The construction of our causal variational approach to initial value problems, relied on the enforcement of different types of conditions on the dynamical degrees of freedom, be it in the form of initial data, spatial boundary data or the connecting condition of the forward and backward branches. So far we enforced these conditions implicitly, which is valid in an analytic setting. In preparation for a numerical treatment, we however have to incorporate these conditions explicitly in the system action. This is possible by use of Lagrange multipliers. I.e we amend <ref> by additional sets of multiplier functions $\lambda,\tilde \lambda$ for initial conditions and $\gamma,\tilde\gamma$ for the connecting conditions. The tilde in $\tilde \lambda$ and $\tilde \gamma$ indicates that these Lagrange multiplier functions are associated with constraints involving temporal derivatives, while the constraints enforced by $\lambda$ and $\gamma$ make reference only to the values. To enforce spatial boundary conditions we introduce the Lagrange multiplier functions $\kappa,\tilde \kappa, \xi,\tilde \xi$, where the tilde refers to those implementing the boundary conditions at $\vec{\sigma}^{\rm f}$ and those without at $\vec{\sigma}^{\rm i}$. Whenever we refer to quantities for which explicit Lagrange multipliers have been taken into account, a superscript letter $\rm L$ will be added. For notational simplicity we here refer only to Dirichlet boundary conditions in spatial direction, which can be generalized to Neumann or mixed boundary conditions

\begin{align}
&\nonumber{\cal E}^{\rm L}_{\rm IBVP}=\int d^{(d+1)}\Sigma\; \Big\{ E_{\rm BVP}[X_1,\partial_a X_1,\phi_1, \partial_a\phi_1] -E_{\rm BVP}[X_2,\partial_a X_2,\phi_2, \partial_a\phi_2]\Big\} \\
\nonumber+&\int \prod_{a=1}^{d}d\Sigma_a\Big\{ \lambda_\mu\big( X_1^\mu(\tau^{\rm i},\vec{\sigma}) - X^\mu_{\rm IC}\big) +\lambda_\phi \big( \phi_1(\tau^{\rm i},\vec{\sigma}) - \phi_{\rm IC}\big)\\
\nonumber+&\tilde \lambda_\mu\big( \partial_0 X_1^\mu(\tau^{\rm i},\vec{\sigma}) - \dot X^\mu_{\rm IC}\big)  +\tilde \lambda_\phi \big( \partial_0\phi_1(\tau^{\rm i},\vec{\sigma}) - \dot \phi_{\rm IC}\big)\\
\nonumber+&\gamma_\mu\big( X_1^\mu(\tau^{\rm f},\vec{\sigma}) - X_2^\mu(\tau^{\rm f},\vec{\sigma})\big) +\gamma_\phi \big( \phi_1(\tau^{\rm f},\vec{\sigma}) - \phi_2(\tau^{\rm f},\vec{\sigma}) \big)\\
+&\tilde \gamma_\mu\big( \partial_0 X_2^\mu(\tau^{\rm f},\vec{\sigma}) - \partial_0 X_2^\mu(\tau^{\rm f},\vec{\sigma})\big)  +\tilde \gamma_\phi \big( \partial_0\phi_1(\tau^{\rm f},\vec{\sigma}) - \partial_0\phi_2(\tau^{\rm f},\vec{\sigma})\big)\Big\}\\
\nonumber+&\sum_{j=1}^{d} \int \prod_{\tiny\begin{array}{c}a=0\\a\neq j\end{array}}^{d} d\Sigma_a\Big\{\kappa^j_{\mu} \big( X_1^\mu(\sigma_j^{\rm i}) - X^\mu_{\rm sBCL}(\sigma_j^{\rm i})\big)
\nonumber+\xi^j_{\mu} \big( X_2^\mu(\sigma_j^{\rm i}) - X^\mu_{\rm sBCL}(\sigma_j^{\rm i})\big)\\
\nonumber+&\tilde \kappa^j_{\mu} \big( X_1^\mu(\sigma_j^{\rm f}) - X^\mu_{\rm sBCR}(\sigma_j^{\rm f})\big)
\nonumber+\tilde \xi^j_{\mu} \big( X_2^\mu(\sigma_j^{\rm f}) - X^\mu_{\rm sBCR}(\sigma_j^{\rm f})\big)\\
\nonumber+&\kappa^j_{\phi} \big( \phi_1(\sigma_j^{\rm i}) - \phi_{\rm sBCR}(\sigma_j^{\rm i})\big)
\nonumber+\xi^j_{\phi} \big( \phi_2(\sigma_j^{\rm i}) - \phi_{\rm sBCR}(\sigma_j^{\rm i})\big)\\
+&\tilde \kappa^j_{\phi} \big( \phi_1(\sigma_j^{\rm f}) - \phi_{\rm sBCL}(\sigma_j^{\rm f})\big)
+\tilde \xi^j_{\phi} \big( \phi_2(\sigma_j^{\rm f}) - \phi_{\rm sBCL}(\sigma_j^{\rm f})\big)\Big\},\label{eq:novelactionEIVPl1}\\
=&\int d^{(d+1)}\Sigma\; \Big\{ E^{\rm L}_{\rm BVP}[X_1,\partial_a X_1,\phi_1, \partial_a\phi_1] - E^{\rm L}_{\rm BVP}[X_2,\partial_a X_2,\phi_2, \partial_a\phi_2]\Big\}.\label{eq:novelactionEIVPlabs}
\end{align}

Above in the second and third line, we enforce the initial conditions on the forward branch of $X^\mu$ and the field $\phi$. I.e. we specify their value and derivative at initial temporal parameter $\tau^{\rm i}$. In the fourth and fifth line we connect the two branches at final temporal parameter $\tau^{\rm f}$. The different quantities $\lambda,\tilde\lambda,\gamma,\tilde \gamma$ are functions of all the $d$ spatial parameters $\vec{\sigma}$ but do not depend on the temporal parameter $\tau$.

Enforcement of the spatial boundary conditions constitutes the rest of the terms, where separate Lagrange multiplier functions are applied to fix left and right spatial boundaries, corresponding to spatial parameters $\vec{\sigma}^{\rm i}$ and $\vec{\sigma}^{\rm f}$ respectively. The spatial boundary conditions are enforced on each of the branches separately with different Lagrange multiplier functions. Here the quantities $\kappa^j,\tilde \kappa^j,\xi^j,\tilde\xi^j$ are functions that depend on the temporal parameter $\tau$ and $(d-1)$ spatial parameters, i.e. they depend on all components of $\vec{\sigma}$, except for the $j$-th parameter $\sigma_j$, indicate by the superscript.

As discussed in [44] for the world-line formalism, the Lagrange multiplier contributions can be absorbed into the action integral by amending each of the boundary terms with delta-functions, placing the multiplier functions at the correct initial or final parameter values. This is indicated in <ref> by the redefinition of  $E_{\rm BVP}\to E^{\rm L}_{\rm BVP}$. The presence of the Lagrange multiplier terms affects the expressions one obtains for the equations of motion, and for the Noether charge in (<ref>) and (<ref>) we have to take derivatives of $E^{\rm L}_{\rm BVP}$ instead of the boundary value Lagrangian $E_{\rm BVP}$. This leads to a modified Noether current

\begin{align}
\partial_a(J^{\rm L})^a=  \partial_a\Big[ \frac{\partial  E^{\rm L}_{\rm BVP}}{\partial (\partial_a X^\mu)} \delta X^\mu \Big] =0,\label{eq:NoetherCurmodL}
\end{align}

where the $a$-th component is denoted $(J^{\rm L})^a$.
In turn we find the following modification of the Noether charge density

\begin{align}
(J^{\rm L})^0=  \Big(\frac{\partial E_{\rm BVP}}{\partial (\partial_0 X^\mu)} + \tilde\lambda_\mu\delta(\tau-\tau^{\rm i}) +\tilde \gamma_\mu\delta(\tau-\tau^{\rm f}) \Big) \delta X^\mu. \label{eq:NoetherChrgDenmodL}
\end{align}

Note that only those Lagrange multipliers $\tilde\lambda_\mu,\tilde\gamma_\mu$ appear, which in the action are associated with constraints involving $\tau$ derivatives of the coordinate maps. We will provide a concrete example of the corrections to the conserved charges when discussing wave propagation in $1+1$ dimensions in <ref>.

At this point we have furnished the full continuum variational formalism in $(d+1)$ dimensions. Elevating coordinate maps to dynamical degrees of freedom led us to the action ${\cal E}_{\rm BVP}$ in <ref>, on which we base our genuine IBVP formulation ${\cal E}_{\rm IVP}$ in <ref> and whose underlying requirements in the form of initial, connecting, and boundary conditions have been made explicit in ${\cal E}^{\rm L}_{\rm IBVP}$ in <ref>. It now remains to discretize this variational approach to IBVPs, which is the focus of the next section.

§ DISCRETIZATION OF THE IBVP

What distinguishes our novel variational approach to IBVPs from the conventional formulation is the presence of dynamical coordinate maps. As we will show below, these maps allow us to discretize the system on the level of the underlying parameters $\Sigma$ instead of the space-time coordinates $X$. In turn we obtain a discretized action  which retains all space-time symmetries and thus preserves the continuum Noether charges. The fact that the coordinate maps adapt to the field dynamics in a way that preserves the continuum symmetries provides a mechanism, where the coarseness of the dynamically emerging space-time grid is automatically adjusted, implementing in effect an automatic adaptive mesh refinement (AMR) procedure.

Retention of symmetries and the preservation of the associated Noether charges in the discrete setting is the central goal of our discretization strategy. The derivation of Noether currents and charges requires the application of integration-by-parts (IBP), hence, our discretization will be based on summation-by-parts (SBP) finite difference operators [22, 23, 24], which mimic IBP exactly in the discrete setting.

Our notation in the discrete setting is as follows: discrete vectors will be referred to as bold letters, such as $\bm X$, $\bm \phi$ or $\bm f$. Matrix operators that act on these vectors are denoted by capital double-struck letters, e.g. $\mathds{H}$ or $\mathds{D}$. In case we wish to refer to operations affecting functions only along a single parameter (either $\tau$ or a $\vec{\sigma}$), we will use instead lowercase double-struck letters, such as $\mathbb{h}_a$ or $\mathbb{d}_a$, where subscripts denote which of the $(d+1)$ space-time dimension the matrix operator is associated with. The identity matrix in $(d+1)$ dimensions is thus referred to as $\mathds{I}$ whereas its one-dimensional counterpart is $\mathbb{i}$.

The underlying $(d+1)$ dimensional parameter space described by $\Sigma^a=(\tau,\vec{\sigma})^a$ is discretized on a hypercubic grid with $N_a$ points and grid spacing $\Delta_a$, not necessarily equal for each dimension. Once we select the intervals $\Sigma_a\in [\Sigma_a^{\rm i},\Sigma_a^{\rm f}]$ in which the parameters lie, the grid spacings follow as $\Delta_a=(\Sigma_a^{\rm f}-\Sigma_a^{\rm i})/(N_a-1)$. Any function  $f(\Sigma)$, which depends on the $\Sigma_a$ parameters, is discretized as an array $ \bm f$ with $= \prod_{a=0}^{d} N_a$ entries. We refer to the discretized dynamical degrees of freedom, the coordinate maps and field on the forward and backward branch, as ${\bm X}^\mu_{1,2}$ and ${\bm \phi}_{1,2}$ respectively. The order of the individual entries in these arrays determines the form of tensor products (implemented as Kronecker products) deployed below. We choose to let $\Sigma_0=\tau$ run slowest and $\Sigma_d=\sigma_d$ fastest.
We need to specify which entry of a discrete array ${\bm f}$ refers to a specific position on the abstract parameter grid. Take e.g. the coordinates $(\tau^{\rm i}+\Delta\tau n_0,\sigma_1^{\rm i}+\Delta\sigma_1 n_1,\ldots,\sigma_d^{\rm i}+\Delta\sigma_d n_d)$, which correspond to taking $n_0$ steps with $\Delta\tau$ and $n_i$ steps with $\Delta\sigma_i$.  Correspondingly we introduce an index function $(\tau^{\rm i}+\Delta\tau n_0,\sigma_1^{\rm i}+\Delta\sigma_1 n_1,\ldots,\sigma_d^{\rm i}+\Delta\sigma_d n_d)= (\ldots(n_0 N_1+n_1)N_2 \ldots + n_d)$ that computes the integer position along the discrete function array corresponding to that specific position. The counter variables $n_0\in[0,N_\tau-1]$ and $n_a\in[0,N_a-1]$ take on only integer values. In $(1+1)$ dimensions, the index is computed explicitly as $(n_\tau,n_\sigma)=(n_\tau N_\sigma +n_\sigma)$.

Integration-by-parts rests on the interplay of differentiation and integration. We start by approximating the $(d+1)$ dimensional integral in the action functional. Note that discretization here refers to the abstract parameters. Let us introduce an inner product on discretized functions ${\bm f}$ and ${\bm g}$, which is characterized by a positive definite diagonal matrix $\mathds{H}\in\mathbb{R}^{\texttt{TotVol}\times\texttt{TotVol}}$, so that $\int d^{(d+1)}\Sigma \; f(\Sigma)g(\Sigma) \approx ( {\bm f}, {\bm g}) = {\bm f}^T \mathds{H} {\bm g}$.
Since in our action, the boundaries of the integrals in each direction $\Sigma_a$ are independent from the values of the other integration variables and we use quadrature rules, which are local to each dimension, we can consider the multidimensional integral described by $\mathds{H}$ as consecutive and independent applications of one-dimensional integration. Hence the matrix $\mathds{H}$ can be written as a tensor product of one-dimensional quadrature rules $\mathbb{h}_a\in\mathbb{R}^{N_a\times N_a}$

\begin{align}
\mathds{H}=\mathbb{h}_0\otimes\mathbb{h}_1\otimes\ldots \otimes\mathbb{h}_d,
\end{align}

where $\mathbb{h}_0$ implements integration along $\tau$ and $\mathbb{h}_i$ along the $\sigma_i$ direction.

Similarly, we define a finite difference operator, which implements the SBP property in one dimension and use the tensor product to  construct the corresponding operator in multiple dimensions. To obtain a finite difference operator $\mathbb{d}_a\in\mathbb{R}^{N_a\times N_a}$ of order $r$ that approximates the continuum derivative $d/d\Sigma_a$ and which mimics IBP exactly, it must fulfill the following SBP properties

\begin{align}
&  \mathbb{d}_a{\bm \Sigma}_a^n = n {\bm \Sigma}_a^{n-1} \quad {\rm for\,n \leq r, \; no\, summation\,in\,a\, implied,\;} {\bm \Sigma}_a{\rm \in\mathbb{R}^{N_a}}\label{eq:derivapprx}\\
& \mathbb{d}_a=\mathbb{h}_a^{-1}\mathbb{q}_a, \quad \mathbb{q}_a^{\rm T}+\mathbb{q}_a=\mathbb{e}_a^N-\mathbb{e}_a^0={\rm diag}[-1,0,\ldots,0,1].\label{eq:SBP}
\end{align}

Here ${\bf \Sigma}_a$ denotes a vector with $N_a$ entries, which linearly increase in magnitude, e.g., ${\bf \Sigma}_a=(0,\Delta_a,2\Delta_a,\ldots)$. The matrix $\mathbb{q}_a$ encodes the so-called SBP property and <ref> establishes that $\mathbb{d}_a$ acts as a derivative, using polynomials as the underlying function space (for a definition of SBP operators on general function spaces see [45, 46, 47]). The matrices $\mathbb{e}^0_a$ and $\mathbb{e}^N_a$ project out the values at the appropriate boundaries of the $\Sigma_a$ direction. To be able to consistently act on our arrays ${\bm X}^\mu_{1,2}$ and ${\bm \phi}_{1,2}$ we must define $\mathds{D}_a\in\mathbb{R}^{\texttt{TotVol}\times\texttt{TotVol}}$ via the tensor product with unit matrices $\mathbb{i}$ as

\begin{align}
\mathds{D}_a= \overbracket{\underbracket{\mathbb{i} \otimes \ldots \otimes \mathbb{i}}_{\rm a\, times} \otimes \mathbb{d}_a \otimes \mathbb{i} \ldots \otimes \mathbb{i}}^{\rm d+1\, times}
\end{align}

As concrete realizations of the $\mathbb{h}_a$ and $\mathbb{d}_a$ we will use the two lowest order diagonal norm SBP schemes referred to as  and . While the former exhibits second order accuracy in the interior, the latter is of fourth order. The accuracy on the boundary is reduced by one order for  and by two orders for [For the  operator this is evident as the second order symmetric stencil in the interior in <ref> reduces to a first order forward or backward stencil on the boundary.]. The quadrature rule for the  scheme is the well-known trapezoidal rule

\begin{equation}
\mathbb{h}_a^{[1,2,1]}=\Delta \Sigma_a \left[ \begin{array}{ccccc} 1/2 & & & & \\ &1 & & &\\ & &\ddots && \\ &&&1&\\ &&&&1/2 \end{array} \right],
\quad
\mathbb{d}_a^{[1,2,1]}=
\frac{1}{2 \Delta \Sigma_a}
\left[ \begin{array}{ccccc} -2 &2 & & &\\ -1& 0& 1& &\\ & &\ddots && \\ &&-1&0&1\\ &&&-2&2 \end{array} \right].\label{eq:SBP121}
\end{equation}

The uniquely defined  operator features the lowest order central symmetric stencil in the interior and the forward and backward stencils at the left and right boundary of the parameter interval respectively.
To achieve fourth order accuracy in the interior and second order on the boundary, the corresponding symmetric stencil is modified by a collection of $4\times 4$ entries in the corners

\begin{align}
\nonumber &\mathbb{h}_a^{[2,4,2]}=\Delta \Sigma_a \left[ \begin{array}{ccccccc}
\frac{17}{48} & & & & & & \\
& \frac{59}{48} & & & & & \\
& &\frac{43}{48} & & & & \\
& & &\frac{49}{48} & & & \\
& & & & &1 & \\
& & & & & &\ddots \\
\end{array} \right],\end{align}

\begin{align}
\frac{1}{\Delta \Sigma_a}
\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}& & \\
\end{array} \right].\label{eq:SBP242}
\end{align}

Note that the left and right eigenvectors of the matrices $\mathbb{d}_a$ are not necessarily the same. In fact it has been shown in [15] that while the right eigenvectors contain the physical zero mode, i.e. the constant function, the zero modes among the left eigenvectors are so-called $\pi$-modes (a manifestation of the doubler problem of lattice field theory [48]). When SBP operators are used in the setting of discretized partial differential equations, only their right eigenvectors play a role. Since in that case only the constant function is annihilated by the derivative operator, one calls $\mathds{D}_a$ null-space consistent (see ref. [49]). In the discretization of an action functional, which is quadratic in the derivatives, additional care must be taken. When locating the critical point of the action, the $\pi$-mode will, if not tamed, contaminate the solution [15]. In this sense the matrices $\mathbb{d}_a$ of <ref> are not yet null-space consistent.

In prior publications [15, 38] it has been established that in a one-dimensional setting, the $\pi$-mode can be avoided by including boundary information into $\mathds{D}_a$, taking inspiration from the technique of simultaneous-approximation terms (SAT) [50]. The strategy relies on incorporating penalty terms related to boundary information into a novel derivative $\bar{\mathds{D}}_a$, without unphysical zero modes[Since the functions we apply the derivative operator on are purely real, the Wilson-term regularization [48] often deployed in lattice field theory is not applicable here.]

\begin{align}
&\bar{\mathds{D}}^f_a {\bm f} = \mathds{D}_a {\bm f} + \sigma_0\mathds{H}^{-1}\mathds{E}^0_a({\bm f} - {\bm f}_{\rm bnd})= \mathds{D}_a {\bm f} + \mathds{S}_a({\bm f} - {\bm f}_{\rm bnd}),\label{eq:regSBP}
\end{align}

where the projector matrices $\mathds{E}^0_a$ and $\mathds{S}_a$ are defined as

\begin{align}
\mathds{E}^0_a= \overbracket{\underbracket{\mathbb{i}\otimes \ldots \otimes \mathbb{i}}_{\rm a\, times} \otimes \mathbb{e}^0_a \otimes \mathbb{i} \ldots \otimes \mathbb{i}}^{\rm d+1\, times},\quad  \mathds{S}_a=\overbracket{\underbracket{\mathbb{i}\otimes \ldots \otimes \mathbb{i}}_{\rm a\, times} \otimes \sigma_0\mathbb{h}_a^{-1}\mathbb{e}^0_a \otimes \mathbb{i} \ldots \otimes \mathbb{i}}^{\rm d+1\, times}
\end{align}

We choose $\sigma_0=1$ to obtain a regularized multidimensional $\bar{\mathds{D}}^f_a$ with eigenvalues featuring a positive definite real part, in agreement with the one-dimensional setting in [15, 38].

As we only consider Dirichlet boundary conditions in spatial directions, the array ${\bm f}_{\rm bnd}$ contains values of the function ${\bm f}$ only on the boundary itself and may be set to zero otherwise. Note that a different regularized SBP operator ensues for each degree of freedom ($\bar{\mathds{D}}^\mu_a$ to act on ${\bm X}^\mu_{1,2}$ and $\bar{\mathds{D}}^\phi_a$ to act on ${\bm \phi}_{1,2}$), as each may be assigned different boundary behavior. The above definition entails that the derivative operator in the temporal direction $\tau$ will include the initial conditions of the corresponding degree of freedom, while for the spatial derivative operator we will use the spatial boundary conditions at $\sigma^{\rm i}$.

The shift operation involving ${\bm f}_{\rm bnd}$ in <ref> can be efficiently implemented by introducing affine coordinates. The introduction of affine coordinates entails that all discrete function arrays are amended by one more entry of value one and each operator matrix is endowed with one more row and column with the value one placed in the lower right corner. The new column on the right contains the shift. The new row, except for the corner is filled with the values zero. The resulting structure of the regularized SBP operator is sketched below with the conventional structures shaded in purple and the affine coordinate additions in orange

\begin{equation}
\bar{\mathds{D}}_a^f \bar{{\bf f}} =\begin{tikzpicture}[baseline={-0.5ex},mymatrixenv]
\matrix [mymatrix,inner sep=4pt] (m)
&    &   &\\
& \mathds{D}_a+\mathds{S}_a & & -\mathds{S}_a{\bf f}_{\rm bnd}   \\
&    &   &\\
& 0   &   & 1\\
\begin{scope}[on background layer,rounded corners]
\node [fit=(m-1-1) (m-3-3),purpleish,inner xsep=1.5pt,inner ysep=2.5pt]{};
\node [fit=(m-1-4) (m-3-4),orangeish,inner xsep=1.5pt,inner ysep=2.5pt]{};
\node [fit=(m-4-1) (m-4-3),orangeish,inner xsep=1.5pt,inner ysep=2.5pt]{};
\node [fit=(m-4-4) (m-4-4),orangeish,inner xsep=1.5pt,inner ysep=2.5pt]{};
\end{scope}
\end{tikzpicture}
\begin{tikzpicture}[baseline={-0.5ex},mymatrixenv]
\matrix [mymatrix,inner sep=4pt] (m)
\\
{\bf f}\\
\\
1 \\
\begin{scope}[on background layer,rounded corners]
\node [fit=(m-1-1) (m-3-1),purpleish,inner xsep=2.5pt,inner ysep=2.5pt]{};
\node [fit=(m-4-1) (m-4-1),orangeish,inner xsep=1.5pt,inner ysep=2.5pt]{};
\end{scope}
\end{tikzpicture}.\label{eq:sketchregSBP}
\end{equation}

For a more comprehensive discussion of the implementation of the affine coordinate regularization we refer the reader to [15, 38].

With appropriately regularized multidimensional SBP operators in place, we proceed to discretize first the novel action ${\cal E}_{\rm BVP}$, as its Lagrangian $E_{\rm BVP}$ constitutes the building block of the genuine causal IBVP formulation. In addition it allows us to highlight the central property of our approach, the retention of continuum space-time symmetries after discretization. Replacing derivatives by regularized SBP operators we obtain

\begin{align}
{\bm g}_{ab}=G_{\mu\nu} (\bar{\mathds{D}}^\mu_a {\bm X}^\mu)\circ (\bar{\mathds{D}}^\nu_b {\bm X}^\nu), \quad {\rm det}[{\bm g}]=\sum_{i_0,\ldots,i_{d}}\epsilon_{i_0\cdots i_{d}} {\bm g}_{0,i_0}\circ\cdots\circ{\bm g}_{d,i_{d}}\label{eq:discrg}
\end{align}

where the symbol $\circ$ refers to pointwise multiplication of discrete arrays and $\epsilon_{i_0\cdots i_{d}}$ denotes the $(d+1)$ dimensional totally anti-symmetric tensor. For ease of notation when dealing with integration over products of more than two functions, we introduce the vector quantity ${\bf h}=\mathds{H}{\bf 1}$, where ${\bf 1}$ denotes the vector with all entries being one. ${\bf h}$ contains the entries of the diagonal matrix $\mathds{H}$ allowing us to write the integral of ${\bf f}$ as $({\bf f},{\bf 1})={\bf f}^T{\bf h}$. Applying the discretization described above to the continuum action ${\cal E}_{\rm BVP}$ in <ref> results in the following expression

\begin{align}
\nonumber \mathds{E}_{\rm BVP}&[{\bm X}_1^\mu,\bar{\mathds{D}}^\mu_a {\bm X}_1^\mu,{\bm \phi}_1, \bar{\mathds{D}}^\phi_a{\bm \phi}_1]=\\
&\frac{1}{2}\Big\{ \Big(\frac{1}{T}V({\bm \phi}_1)-1\Big)\circ{\rm det}[{\bm g}_1] +\frac{1}{T} (\bar{\mathds{D}}^\phi_a{\bm \phi}_1) \circ  (\bar{\mathds{D}}^\phi_b{\bm \phi}_1)\circ{\rm adj}[{\bm g}_1]_{ab}\Big\}^T {\bm h}.\label{eq:discrEBVP}
\end{align}

Let us start by investigating the discretized induced metric ${\bm g}_{ab}$ in (<ref>). We see that the discrete derivatives in (<ref>) are combined with the space-time metric $G$ in the same way as in the continuum. By Einstein summation convention we sum over both $\mu$ and $\nu$ indices. This entails that under a Lorentz transform, a transformation matrix $\Lambda$ is applied to each of the coordinate maps ${\bf X}$. The discretized expression remains invariant by definition of $\Lambda$, i.e $G_{\mu\nu}\Lambda^\mu_\kappa\Lambda^\nu_\gamma=G_{\kappa\gamma}$. Thus each entry of the matrix ${\bm g}$ individually stays invariant under Lorentz transformations in the discrete setting. Furthermore the discrete derivatives $\bar{\mathds{D}}_a^\mu$ annihilate the constant function exactly, such that adding a constant global shift to each entry of the arrays ${\bf X}^\mu$ will not affect the discrete ${\bm g}$, making it invariant under the full Poincaré group of transformations. In turn, as the only reference to space-time coordinates in <ref> is made via ${\bm g}$ we can conclude that the discretized action $\mathds{E}_{\rm BVP}$ retains the full manifest Poincaré invariance of the continuum action.

This situation is crucially different from the conventional discretization in the space-time coordinates themselves. In our novel formulation it is the underlying abstract parameters $\Sigma$, which are discretized, leaving the entries of the ${\bf X}^\mu$ arrays to take on arbitrary values. In particular the values of the entries of ${\bf X}^\mu$ can be changed by arbitrarily infinitesimal amounts, the prerequisite for an application of Noether's theorem.

In the derivation of the continuum Noether current only integration by parts was deployed and our discretization via SBP operators mimics IBP exactly in the discrete setting. We may thus directly infer that the expression for the Noether current ${\bm J}^a$ and its associated Noether charge density ${\bm q}\equiv{\bm J}^0$ in the discrete setting reads

\begin{align}
\mathds{D}_a{\bf J}^a=  \mathds{D}_a\Big[ \frac{\partial \mathds{E}_{\rm BVP}}{\partial (\mathds{D}_a {\bf X}^\mu)} \delta {\bf X}^\mu \Big] =0, \quad {\bm q} = \frac{\partial \mathds{E}_{\rm BVP}}{\partial (\mathds{D}_0 {\bf X}^\mu)}\delta {\bf X}^\mu.\label{eq:NoetherCurDiscr}
\end{align}

Note that the conserved Noether charge ${\bm Q}$ is obtained from integration of the associated charge density ${\bm q}$ over all spatial direction (c.f. <ref>).

The differentiation with respect to the array containing the derivative of the space-time coordinate maps $\partial/\partial (\mathds{D}_a {\bf X}^\mu)$ is to be understood in a point-wise fashion. In the discrete setting it can be explicitly implemented using the discrete Dirac delta function $ \mathfrak{d}_k$ (also known as lifting operator [51]). When ${\bf f}$ is integrated over $ \mathfrak{d}[k]$ the result must correspond to the $k$-th entry of ${\bf f}$, i.e. ${\bf f}^{\rm T}{\mathds{H}}\mathfrak{d}[k] ={\bf f}[k]$. We can consider both the discrete counterpart to the delta function collapsing integration in all dimensions $\delta^{(d+1)}(\Sigma-\Sigma_{\rm ref})$ and those acting only in a single dimension $\delta(\Sigma^a-\Sigma^a_{\rm ref})$

\begin{align}
\mathfrak{d}[k] = {\mathds{H}}^{-1} {\bf e}_k, \quad \mathfrak{d}^a[j] =  \mathbb{h}_a^{-1} {\bf e}^a_j \label{eq:defdeltas}
\end{align}

where ${\bf e}\in\mathbb{R}^{\texttt{TotVol}}$ and ${\bf e}^a\in\mathbb{R}^{N_a}$ refer to appropriate unit vectors filled with zero except for the value one at the position $k$ or $j$ at which the delta function is non-zero.

Using the discrete delta function defined above, we can mimic continuum functional derivatives $\frac{\delta }{\delta f(z)}\int dx f(x) g(x) = \int dx \delta(x-z)g(x) = g(z)$ in the discrete setting as  $\partial ({\bf f},{\bf g})/\partial {\bf f}[k]= (\mathfrak{d}[k],{\bf g})={\bf g}[k]$.

In preparation for expressing the Noether charge in the discrete setting, we also need to define a matrix operator, which implements spatial integration. I.e. it must transform a discrete array ${\bf f}$ with  entries to an array that only contains $N_\tau$ entries, an operation represented via a $N_\tau\times$ matrix we call $\mathds{H}_\sigma$. For convenient implementation in practice, we define the $\prod_{a>0}^{d} N_a$ component vector $({\bm h}_\sigma)_i=(\mathbb{h}_1\otimes\ldots\otimes\mathbb{h}_d)_{ii}$, which contains the diagonal entries of the spatial quadrature matrix $\mathbb{h}_1\otimes\ldots\otimes\mathbb{h}_d$. Since we decided that our coordinates are ordered in the discrete setting such that spatial coordinates run fastest, the matrix $\mathds{H}_\sigma$ can be written as

\begin{align}
\mathds{H}_\sigma = \left(\begin{array}{cccc} {\bm h}^T_\sigma & 0 & \cdots &0\\
0 & {\bm h}^T_\sigma& \cdots &0\\
\vdots&&&\\
0 & 0 &\cdots& {\bm h}^T_\sigma\end{array}\right).\label{eq:spatint}
\end{align}

In <ref> we showed that in order to explicitly include the necessary conditions related to initial and boundary data, as well as the connection of the forward and backward branch, we can use Lagrange multiplier functions absorbed into a redefinition of the BVP functional $E_{\rm BVP}\to E^{\rm L}_{\rm BVP}$. In the discrete setting we similarly incorporate these terms in a new $\mathds{E}_{\rm BVP}\to {\mathds{E}}^{\rm L}_{\rm BVP}$ using the discrete delta function defined above. The Lagrange multiplier terms affect the resulting equations of motion and the Noether charge. As discussed in detail for the one-dimensional case in [44], the modification to the Noether charge can be uniquely related to the Lagrange multiplier terms associated with the derivatives of the spatial coordinates. Denoting the discretized Lagrange multiplier functions by bold-face Greek letters, the discretized version of the continuum Noether charge density of <ref> reads

\begin{align}
{\bm q}^{\rm L} = \frac{\partial {\mathds{E}}^{\rm L}_{\rm BVP}}{\partial (\mathds{D}_0 {\bf X}^\mu)}\delta {\bf X}^\mu = \Big( \frac{\partial {\mathds{E}}_{\rm BVP}}{\partial (\mathds{D}_0 {\bf X}^\mu)} + \tilde{\bm \lambda}_\mu \circ \mathfrak{d}^0[0] + \tilde{\bm \gamma}_\mu\circ \mathfrak{d}^0[N_0] \Big)\delta {\bf X}^\mu,\label{eq:NoetherCurDiscrLagr}
\end{align}

while the following explicit expression for the corresponding Noether charge ensues

\begin{align}
{\bm Q}^{\rm L} =  \Big( \mathds{H}_\sigma \frac{\partial {\mathds{E}}_{\rm BVP}}{\partial (\mathds{D}_0 {\bf X}^\mu)} + ({\bm h}_\sigma^T \tilde{\bm \lambda}_\mu) \mathfrak{d}^0[0] + ( {\bm h}_\sigma^T \tilde{\bm \gamma}_\mu) \mathfrak{d}^0[N_0] \Big)\delta {\bf X}^\mu.\label{eq:NoetherChargeDiscrLagr}
\end{align}

When considering the concrete implementation of our discretization strategy for wave propagation in $(1+1)$ dimensions in the following section, we will show that the Noether charge density ${\bm q}^{\rm L}$ indeed leads to a Noether charge ${\bm Q}^{\rm L}$ which is exactly conserved over the whole duration of the simulation.

The Lagrange multiplier terms also affect the form of the equations of motion, leading to deviations from the continuum form at the last and second to last $\tau$ slice. However, we will show that, as in the world-line formalism, these modifications do not spoil the convergence of the solution to the correct continuum limit nor do they interfere with the conservation of the Noether charge.

§ PROOF-OF-PRINCIPLE: WAVE PROPAGATION IN $(1+1)$ DIMENSIONS

Having established the general variational approach to IBVPs with dynamical coordinate maps in $(d+1)$ dimensions in the previous sections, we proceed to present an explicit implementation for a dynamical system of central interest, scalar wave propagation in $(1+1)$ dimensional flat space-time. It serves as a demonstration for how our novel approach addresses two and connects to the third of the three central challenges to discretized IBVPs. In the context of the first challenge, we show how space-time symmetries can be retained after discretization and that the continuum Noether charges stay exactly preserved in the discrete context. Taking aim at the second challenge, we show how a non-constant discretization of the space-time coordinates emerges dynamically, realizing a form of automatic adaptive mesh refinement (AMR). Along the way we will see how the presence of the coordinate maps offers new freedom in the implementation of boundary conditions, connecting to the third challenge.

In order to keep computational cost for this proof-of-principle to a minimum, we consider $(1+1)$ dimensions, where the abstract parameter space is two dimensional $\Sigma^a=(\tau,\sigma)^a$ and the most general coordinate maps read $X^0(\tau,\sigma)=t(\tau,\sigma)$ and $X^1(\tau,\sigma)=x(\tau,\sigma)$. In $(1+1)$ dimensions the continuum theory features three space-time symmetries: translations in time $t$, translations in space $x$, and relativistic rotations involving both the time and space coordinate, so called boosts. I.e. even though spatial rotation are absent in $(1+1)$ dimensions, our proof-of-principle contains a space-time symmetry structure that is similar to the complexity found in higher dimensions. The value of the scale $T$ is taken to be large compared to the energy density stored in the initial conditions of the field degrees of freedom.

In order to avoid complications in numerically locating the critical point of the IBVP action,
we will consider a trivial mapping in space $x(\tau,\sigma)=\sigma$ while leaving the temporal mapping fully dynamic $t=t(\tau,\sigma)$. This choice allows us to focus on the preservation of time translation symmetry after discretization, associated with the concept of a conserved energy.
Conservation of energy, while not only a central physical tenet of the classical dynamics, is vital to the stability of the discretization.

§.§ Continuum formulation

In flat (Minkowski-)space-time in $(1+1)$ dimensions, the metric tensor takes the form $G={\rm diag}[c^2,-1]$. In the presence of dynamical coordinate maps we must determine the corresponding induced metric $g$ on the abstract parameter manifold spanned by the $(\tau,\sigma)$ directions. Let us denote derivatives with respect to $\tau$ as dots, as in $\dot t\equiv\partial t/\partial\tau $ and $\dot x\equiv \partial x/\partial\tau $, and we will indicate derivatives w.r.t. the spatial parameter $\sigma$ as primes $t^\prime\equiv \partial t/\partial\sigma$ and $ x^\prime\equiv \partial x/\partial\sigma$. The explicit expression for the Jacobian of the coordinate maps then is

\begin{align}
J^\mu_a=\partial X^\mu/\partial \Sigma^a=\left(\begin{array}{cc} \dot t& t^\prime \\ \dot x& x^\prime \end{array}\right).
\end{align}

The induced metric $g=J^TGJ$ and its adjugate thus read

\begin{align}
g=\left( \begin{array}{cc} c^2\td^2-\xd^2 & c^2\td\tp - \xd\xp \\c^2 \tp\td - \xp\xd & c^2(\tp)^2-(\xp)^2\end{array}\right), \; {\rm adj}[g]= \left( \begin{array}{cc} c^2(\tp)^2-(\xp)^2 & \xd\xp-c^2\td\tp\\ \xp\xd-c^2 \tp\td & c^2\td^2-\xd^2\end{array}\right),
\end{align}

and the determinant reduces to

\begin{align}
{\rm det}[g]= -c^2(\td\xp-\xd\tp)^2\label{eq:fulldetg}.
\end{align}

It is remarkable that the full determinant reduces to such a simple expression.

Setting the potential term to zero as $V(\phi)=0$ in our novel action ${\cal S}_{\rm BVP}$ (introduced in <ref>) corresponds to a system that describes the propagation of scalar waves with velocity $c$. Derivatives of the dynamical field $\phi(\tau,\sigma)$ are denoted similarly to those of the coordinate maps as $\dot \phi\equiv\partial \phi/\partial\tau $ and $\phi^\prime\equiv \partial \phi/\partial\sigma$. For the $(1+1)$ dimensional setting we thus have

\begin{align}
\nonumber {\cal S}_{\rm BVP}=&\int d\tau d\sigma \, \big(-T\big) \sqrt{ \Big(-{\rm det}[g] + \frac{1}{T} \partial_a\phi(\Sigma)\partial_b\phi(\Sigma) {\rm adj}[g]_{ab}}\\
=&\int d\tau d\sigma \, \big(-T\big) \Big\{ c^2(\td\xp-\xd\tp)^2\label{eq:novelaction1p1}\\
\nonumber+&\frac{1}{T} \Big( \pd^2(c^2(\tp)^2-(\xp)^2) + 2\pd\pp(\xd\xp-c^2\td\tp) + (\pp)^2(c^2\td^2-\xd^2) \Big) \Big\}^{1/2}.
\end{align}

Note how each term under the square root contains exactly two $\tau$ and two $\sigma$ derivatives. This means that the integral is not only invariant under a differentiable redefinition of $(\tau,\sigma)\to f(\tau,\sigma)$ but that at the same time the dimensions of the parameters $(\tau,\sigma)$ are irrelevant; the dimensions of the parameters cancel between the measure and the derivative terms. Whenever temporal and spatial coordinate maps appear in <ref> they do so quadratically with a relative minus sign, which encodes the proper invariance of the expressions under Lorentz transformations. As only dotted and primed coordinates contribute, global shifts in time and space too constitute manifest symmetries.
In the following we will, for notational clarity, choose natural units such that $c=1$. Specifying to the trivial spatial mapping $x(\tau,\sigma)=\sigma$ and removing the square root from the integrand we obtain the following action on which we base our proof-of-principle

\begin{align}
{\cal E}_{\rm BVP}\overset{x=\sigma}{=}\int &d\tau d\sigma \, \frac{1}{2} \Big\{ (\td)^2+\frac{1}{T} \Big( \pd^2((\tp)^2-1) - 2 \pd\pp\td\tp + (\pp)^2(\td^2) \Big) \Big\}. \label{eq:novelactionE1p1}
\end{align}

Note that even with the trivial mapping $x=\sigma$ we still retain manifest symmetry under global time translation, since each reference to the $t$-mapping in <ref> involves either its $\tau$ or $\sigma$ derivative. Let us derive the equations of motion for this system according to (<ref>).

The variation of this action leads to

\begin{align}
\nonumber\delta {\cal E}_{\rm BVP}= \int d\tau d\sigma \Big\{& \td\delta\td + \frac{1}{T} \Big( \pd\delta\pd((\tp)^2-1) +\pd^2\tp\delta\tp - \delta\pd\pp\td\tp - \pd\delta\pp\td\tp\\
& -\pd\pp\delta\td \tp - \pd\pp\td\delta\tp + \pp\delta\pp\td^2+(\pp)^2\td\delta\td\Big)\Big\}.
\end{align}

Carrying out integration by parts results in

\begin{align}
\delta {\cal E}_{\rm BVP}= &-\int d\tau d\sigma \Big\{\Big[ \ddot t + \frac{1}{T}\Big\{ \frac{\partial}{\partial \sigma}\Big( \pd^2\tp - \pd\pp\td \Big) + \frac{\partial}{\partial \tau} \Big( (\pp)^2\td - \pd\pp\tp \Big)\Big\}\Big] \delta t\label{eq:1p1dvareomt}\\
&+ \Big[ \frac{\partial}{\partial \tau} \Big( \pd ((\tp)^2-1) -\pp\td\tp \Big) + \frac{\partial}{\partial \sigma} \Big( \pp \td^2 - \pd\td\tp \Big) \Big]\delta \phi \Big\}\label{eq:1p1dvareomphi}\\
\nonumber &+\int d\sigma \Big\{ \Big[ \dot t + \Big( (\pp)^2\td - \pd\pp\tp \Big)\Big] \delta t\Big\} \Big|_{\tau^{\rm i}}^{\tau^{\rm f}} + \int d\tau \Big\{ \Big[ \pd^2\tp - \pd\pp\td\Big]\delta t   \Big\}\Big|_{\sigma^{\rm i}}^{\sigma^{\rm f}}\\
\nonumber &+\int d\sigma \Big\{ \Big[ \pd ((\tp)^2-1) -\pp\td\tp \Big] \delta \phi\Big\} \Big|_{\tau^{\rm i}}^{\tau^{\rm f}} + \int d\tau \Big\{ \Big[ \pp \td^2 - \pd\td\tp  \Big]\delta \phi   \Big\}\Big|_{\sigma^{\rm i}}^{\sigma^{\rm f}},
\end{align}

from which we can read off the coupled equations of motion for $t$ and $\phi$ from (<ref>) and (<ref>) respectively

\begin{align}
\ddot t + \frac{1}{T}\Big\{ \frac{\partial}{\partial \sigma}\Big( \pd^2\tp - \pd\pp\td \Big) + \frac{\partial}{\partial \tau} \Big( (\pp)^2\td - \pd\pp\tp \Big)\Big\}&=0,\label{eq:1p1deomt}\\
\frac{\partial}{\partial \tau} \Big( \pd ((\tp)^2-1) -\pp\td\tp \Big) + \frac{\partial}{\partial \sigma} \Big( \pp \td^2 - \pd\td\tp \Big)&=0.\label{eq:1p1deomphi}
\end{align}

We may recover the conventional wave equation from <ref> after re-instituting the explicit factors of $c$, by restricting to a trivial time map $t(\tau,\sigma)=\tau$

\begin{align}
0&\quad=\quad \frac{\partial}{\partial \tau} \Big( \pd (c^2(\tp)^2-1) -c^2\pp\td\tp \Big) + \frac{\partial}{\partial \sigma} \Big( c^2 \pp \td^2 - c^2\pd\td\tp \Big)\\
&\overset{t=\tau, x=\sigma}{=}-\frac{\partial}{\partial t} \Big( \pd \Big) + \frac{\partial}{\partial x} \Big( c^2 \pp \Big).
\end{align}

In the limit of $T\to\infty$ we would recover from <ref> the simple relation $\ddot t(\tau,\sigma)=0$ for the time mapping, which is solved by a sheet, constant in the $\sigma$ direction and with constant slope $\dot t_{\rm IC}$ along $\tau$. Note that for $T\to\infty$ the action only contains reference to $\td^2$, meaning that it decomposes into independent temporal maps for each value of $\sigma$. As there are no terms involving spatial derivatives $\tp$ in the action, also no spatial boundary terms arise when deriving the governing equations. In the presence of fields the dynamics of $\phi$ couples to both $\tau$ and $\sigma$ derivatives of $t$, leading to non-trivial dynamics of the time mapping.

In contrast to the conventional wave equation, the equation of motion (<ref>) for the field here not only contains second derivatives in $\tau$ and $\sigma$ but due to the dynamical nature of the $t$ mapping, also mixed terms.

Let us take a look at the spatial boundary conditions necessary for achieving equivalence between the equations of motion <ref> and the critical point of the action <ref>. Note that the temporal boundary conditions will be taken care of by the forward backward construction, when going over from ${\cal E}_{\rm BVP}\to{\cal E}_{\rm IBVP}$.
We thus find that the following two terms must vanish

\begin{align}
\int d\tau \Big\{ \Big[ \pd^2\tp - \pd\pp\td\Big]\delta t   \Big\}\Big|_{\sigma^{\rm i}}^{\sigma^{\rm f}}\overset{!}{=}0, \qquad \int d\tau \Big\{ \Big[ \pp \td^2 - \pd\td\tp  \Big]\delta \phi   \Big\}\Big|_{\sigma^{\rm i}}^{\sigma^{\rm f}}\overset{!}{=}0.\label{eq:bndcnd1p1}
\end{align}

In contrast to the conventional formulation of IBVPs we now have both $\phi$ and $t$ appearing in these boundary expressions. This amounts to a new freedom available to us in order to construct boundary conditions, such that the boundary terms vanish.

In our proof-of-principle here, we simply choose vanishing spatial Dirichlet boundary conditions for the field $\phi$. This entails that $\phi$ and $\delta \phi$, as well as $\dot \phi$ (due to the time independence of Dirichlet boundary conditions), are all zero on the spatial boundary. This choice leaves us with vanishing boundary terms in <ref> without the need to specify a boundary condition on $t$. Conversely one could imagine boundary conditions enforced for $t$, which would make (<ref>) vanish without the need to provide an explicit boundary condition on $\phi$.
I.e. the presence of the coordinate maps provides a unique and novel flexibility in the treatment of boundary conditions for the dynamical fields, especially interesting for non-reflecting boundary conditions.

Having confirmed that for our choice of boundary conditions (<ref>) vanishes, we have established the equations of  motion for our $(1+1)$ dimensional wave propagation.

Following <ref>, the components of the conserved Noether current $\partial_a J_t^a = 0$ associated with manifest invariance under infinitesimal time translations $\delta t$, i.e. $\delta X^\mu=\delta t \delta^{\mu 0}$ read

\begin{align}
&J_t^0=\frac{\partial E_{\rm BVP}}{\partial (\td) } = \td + \frac{1}{T}\Big( (\pp)^2\td - \pd\pp\tp \Big)\\
&J_t^1=\frac{\partial E_{\rm BVP}}{\partial (\tp)} =  \frac{1}{T}\Big( (\pd)^2\tp - \pd\pp\td \Big)
\end{align}

The zeroth component of a conserved current defines the Noether charge according to <ref>, which yields

\begin{align}
Q_t(\tau)=\int d\sigma J_t^0 = \int d\sigma \Big\{ \td + \frac{1}{T}\Big( (\pp)^2\td - \pd\pp\tp \Big)\Big\}.\label{eq:1p1dnoethcont}
\end{align}

This quantity in the presence of dynamical coordinate maps constitutes the energy, which depends both on contributions from the fields and from the coordinate maps themselves.

Let us confirm that our choice of boundary condition conserves the Noether charge

\begin{align}
\nonumber \partial_\tau Q_t &\overset{{\rm def.}\, \cref{eq:1p1dnoethcont}}{=} \int d\sigma \frac{\partial }{\partial \tau} \Big\{ \td + \frac{1}{T}\Big( (\pp)^2\td - \pd\pp\tp \Big)\Big\}\\
\nonumber &\overset{ \partial_0J_t^0=-\partial_1J_t^1}{=} -\int d\sigma \frac{\partial }{\partial \sigma} \Big\{ \frac{1}{T}\Big( \pd^2\tp - \pd\pp\td\Big) \Big\}\\
&\qquad = -\frac{1}{T}\Big( \pd^2\tp - \pd\pp\td\Big)\Big|_{\sigma_{\rm i}}^{\sigma^{\rm f}}=0.
\end{align}

In the second line we have exploited that the Noether current is conserved by construction. And indeed due to the vanishing of $\phi$ on the spatial boundaries and thus $\dot \phi$ being zero there, the Noether charge does not change over time.

Note that the Noether charge of <ref>, in contrast to its counterpart in the world-line formalism, makes reference to all dynamical degrees of freedom, i.e. both their spatial and temporal derivatives. This places tight constraints on the growth of derivatives, and while not a proof, suggests that the evolution of the system is stable. This will be considered in future work.

When we wish to formulate a causal IBVP we have to go over to the forward backward construction and incorporate initial and boundary conditions, as well as the connecting conditions for both the field and temporal mapping via explicit Lagrange multipliers. Note that the evolution of the time mapping too is handled as an IBVP, meaning we do not specify an initial and final time but instead an initial time $t_{\rm IC}$ and the derivative of time $\dot t_{\rm IC}$ with respect to the temporal parameter $\tau$. Similarly we have to supply $\phi_{\rm IC}$ and $\dot \phi_{\rm IC}$ at the initial $\tau$ slice. At this point we also add all possible Lagrange multiplier functions $\kappa,\tilde \kappa,\xi, \tilde \xi$ related to the spatial (Dirichlet) boundary conditions

\begin{align}
{\cal E}^{\rm L}_{\rm IBVP}=&\int d\tau d\sigma \, \frac{1}{2} \Big\{ (\td_1)^2+\frac{1}{T} \Big( \pd^2_1((\tp_1)^2-1) - 2 \pd_1\pp_1\td_1\tp_1 + (\pp_1)^2(\td_1^2) \Big) \Big\}\label{eq:1p1dEIBVPc}\\
\nonumber -&\int d\tau d\sigma \, \frac{1}{2} \Big\{ (\td_2)^2+\frac{1}{T} \Big( \pd_2^2((\tp_2)^2-1) - 2 \pd_2\pp_2\td_2\tp_2 + (\pp_2)^2(\td_2^2) \Big) \Big\}\\
\nonumber +&\int d\sigma \Big\{ \lambda^t(\sigma)\Big(t_1(\tau^{\rm i},\sigma)-t_{\rm IC}(\sigma)\Big)+ \lambda^\phi(\sigma)\Big(\phi_1(\tau^{\rm i},\sigma)-\phi_{\rm IC}(\sigma)\Big)\Big\}\\
\nonumber +&\int d\sigma \Big\{ \tilde \lambda^t(\sigma)\Big(\dot t_1(\tau^{\rm i},\sigma)-\dot t_{\rm IC}(\sigma)\Big)+ \tilde \lambda^\phi(\sigma)\Big(\dot \phi_1(\tau^{\rm i},\sigma)-\dot \phi_{\rm IC}(\sigma)\Big)\Big\}\\
\nonumber +&\int d\sigma \Big\{ \gamma^t(\sigma)\Big(t_1(\tau^{\rm f},\sigma)-t_2(\tau^{\rm f},\sigma)\Big)+ \gamma^\phi(\sigma)\Big(\phi_1(\tau^{\rm f},\sigma)-\phi_2(\tau^{\rm f},\sigma)\Big)\Big\}\\
\nonumber +&\int d\sigma \Big\{ \tilde\gamma^t(\sigma)\Big(\dot t_1(\tau^{\rm f},\sigma)-\dot t_2(\tau^{\rm f},\sigma)\Big)+ \tilde \gamma^\phi(\sigma)\Big(\dot \phi_1(\tau^{\rm f},\sigma)-\dot \phi_2(\tau^{\rm f},\sigma)\Big)\Big\}\\
\nonumber +&\int d\tau \Big\{ \kappa^{t}(\tau)\Big(t_1(\tau,\sigma^{\rm i})-t_{\rm sBCL}(\tau)\Big)+ \kappa^{\phi}(\tau)\Big(\phi_1(\tau,\sigma^{\rm i})-\phi_{\rm sBCL}(\tau)\Big)\Big\}\\
\nonumber +&\int d\tau \Big\{ \tilde\kappa^t(\tau)\Big(t_1(\tau,\sigma^{\rm f})-t_{\rm sBCR}(\tau)\Big)+ \tilde\kappa^\phi(\tau)\Big(\phi_1(\tau,\sigma^{\rm f})-\phi_{\rm sBCR}(\tau)\Big)\Big\}\\
\nonumber +&\int d\tau \Big\{ \xi^{t}(\tau)\Big(t_2(\tau,\sigma^{\rm i})-t_{\rm sBCL}(\tau)\Big)+ \xi^{\phi}(\tau)\Big(\phi_2(\tau,\sigma^{\rm i})-\phi_{\rm sBCL}(\tau)\Big)\Big\}\\
\nonumber +&\int d\tau \Big\{ \tilde\xi^t(\tau)\Big(t_2(\tau,\sigma^{\rm f})-t_{\rm sBCR}(\tau)\Big)+ \tilde\xi^\phi(\tau)\Big(\phi_2(\tau,\sigma^{\rm f})-\phi_{\rm sBCR}(\tau)\Big)\Big\}.
\end{align}

The presence of Lagrange multipliers leads to a modification of the expression for the energy. In the physical limit the values of the degrees of freedom on the forward and backward branch must agree, hence we define the Noether charge simply from the forward time mapping and field

\begin{align}
Q_t^{\rm L}(\tau)=\int d\sigma \Big\{ \td_1 + \frac{1}{T}\Big( (\pp_1)^2\td_1 - \pd_1\pp_1\tp_1 \Big) + \tilde\lambda^t\delta(\tau-\tau^{\rm i})  + \tilde\gamma^t\delta(\tau-\tau^{\rm f}) \Big\}.\label{eq:1p1dnoethcontl}
\end{align}

With the full continuum IBVP action defined in <ref> we are now ready to discretize the action and obtain the classical trajectory by numerically determining its critical point.

§.§  Discretized formulation

In the following we discretize the two-dimensional parameter space $(\tau,\sigma)$ using a grid with $(N_\tau,N_\sigma)$ points. We choose finite intervals for the values of $\tau\in[\tau^{\rm i},\tau^{\rm f}]$ and $\sigma\in[\sigma^{\rm i},\sigma^{\rm f}]$ leading to $\Delta \tau= (\tau^{\rm f}-\tau^{\rm i})/(N_\tau-1)$ and $\Delta \sigma= (\sigma^{\rm f}-\sigma^{\rm i})/(N_\sigma-1)$. The arrays to house the dynamical degrees of freedom are thus ${\bf t}_{1,2}$ and ${\bm \phi}_{1,2}$ each with $=N_\tau N_\sigma$ entries[The spatial mapping was chosen to be trivial $x=\sigma$ and thus does not require to be treated as a separate array.]. For reference to subsets of these arrays let us introduce projection operators $\mathds{P}_a^k$. When applied to an array of  entries, they produce the $k$-th slice in the $a$ direction, e.g. $\mathds{P}_\tau^0[{\bm \phi}]=(\phi(0,0), \phi(0,\Delta\sigma),\ldots,\phi(0,\sigma^{\rm f}))\in \mathbb R^{N_\sigma}$.

As discussed in <ref>, we deploy summation by parts operators $\mathds{D}_\tau=\mathbb{d}_\tau\otimes{\mathbb{i}}$ and $\mathds{D}_\sigma={\mathbb{i}}\otimes\mathbb{d}_\sigma$, as well as approximate integration by the compatible quadrature matrix $\mathds{H}=\mathbb{h}_\tau\otimes\mathbb{h}_\sigma$. The array with the diagonal entries of $\mathds{H}$ we refer to as ${\bm h}$, which is used to express integration of a function $\bm f$ as $({\bm f},{\bm 1})={\bm f}^T \mathds{H} {\bm 1}= {\bm f}^T {\bm h}$. For the temporal direction, where we have a genuine initial value problem, as the solution at the final point $\tau^{\rm f}$ is a priori unknown, we introduce the previously discussed regularization in the SBP operators. This leads to $\bar{\mathds{D}}_{\tau}^t$ and $\bar{\mathds{D}}_{\tau}^\phi$, whose regularization is based on penalty terms including reference to the initial values. We found that when there are spatial Dirichlet boundary conditions enforced, we do not need to introduce regularization, as the $\pi$-mode is avoided automatically by requiring the function to take a specific value at both ends of the interval.

The Lagrange multiplier arrays are referred to as ${\bm \lambda}^t,\tilde{{\bm \lambda}}^t,{\bm \lambda}^\phi,\tilde{{\bm \lambda}}^\phi$ for the initial conditions stored in ${\bm t}_{\rm IC},{\bm \phi}_{\rm IC}$ and $\dot {\bm t}_{\rm IC},\dot {\bm \phi}_{\rm IC}$ arrays. Note that these arrays only contain $N_\sigma$ entries, as they make reference solely to the first $\tau$ slice. Similarly we have ${\bm \gamma}^t,\tilde{{\bm \gamma}}^t,{\bm \gamma}^\phi,\tilde{{\bm \gamma}}^\phi$ for the connecting conditions. As we saw in <ref> if we supply vanishing spatial Dirichlet boundary data to the field on the forward and backward branch using $N_\tau$ component arrays ${\bm \kappa}^\phi,\tilde{{\bm \kappa}}^\phi,{\bm \xi}^\phi,\tilde{{\bm \xi}}^\phi$ no spatial boundary conditions need to be supplied for ${\bm t}$.

We arrive at the following discretized action functional

\begin{align}
\nonumber \mathds{E}_{\rm IBVP}^{\rm L} =& \frac{1}{2}\Big\{ (\bar{\mathds{D}}_\tau^t {\bm t}_1)^2+\frac{1}{T} \Big( (\bar{\mathds{D}}_\tau^\phi {\bm \phi}_1)^2\circ\big( (\mathds{D}_\sigma {\bm t}_1)^2 - 1 \big)\\
\nonumber &\qquad -2(\bar{\mathds{D}}_\sigma^\phi {\bm \phi}_1)\circ(\bar{\mathds{D}}_\tau^\phi {\bm \phi}_1)\circ (\bar{\mathds{D}}_\tau^t {\bm t}_1)\circ(\mathds{D}_\sigma^t {\bm t}_1)  + (\bar{\mathds{D}}_\sigma^\phi {\bm \phi}_1)^2\circ (\bar{\mathds{D}}_\tau^t {\bm t}_1)^2 \Big)\Big\}^T {\bm h}\\
\nonumber-&\frac{1}{2}\Big\{ (\bar{\mathds{D}}_\tau^t {\bm t}_2)^2+\frac{1}{T} \Big( (\bar{\mathds{D}}_\tau^\phi {\bm \phi}_2)^2\circ \big( (\mathds{D}_\sigma {\bm t}_2)^2 - 1 \big)\\
\nonumber &\qquad -2(\bar{\mathds{D}}_\sigma^\phi {\bm \phi}_2)\circ(\bar{\mathds{D}}_\tau^\phi {\bm \phi}_2)\circ (\bar{\mathds{D}}_\tau^t {\bm t}_2)\circ(\mathds{D}_\sigma^t {\bm t}_2) + (\bar{\mathds{D}}_\sigma^\phi {\bm \phi}_2)^2 \circ(\bar{\mathds{D}}_\tau^t {\bm t}_2)^2\Big)\Big\}^T {\bm h}\\
\nonumber+& \big( {\bm \lambda}^t \big)^T \mathbb{h_\sigma} \big( \mathds{P}_\tau^0[{\bm t}_1] - {\bm t}_{\rm IC} \big) + \big( {\bm \lambda}^\phi \big)^T \mathbb{h_\sigma} \big( \mathds{P}_\tau^0[{\bm \phi}_1] - {\bm \phi}_{\rm IC} \big)\\
\nonumber+& \big( \tilde{{\bm \lambda}}^t \big)^T \mathbb{h_\sigma} \big( \mathds{P}_\tau^0[(\mathds{D}_\tau{\bm t}_1)] - \dot{\bm t}_{\rm IC} \big) + \big( \tilde{\bm \lambda}^\phi \big)^T \mathbb{h_\sigma} \big( \mathds{P}_\tau^0[(\mathds{D}_\tau{\bm \phi}_1)] - \dot{\bm \phi}_{\rm IC} \big)\\
\nonumber+& \big( {\bm \gamma}^t \big)^T \mathbb{h_\sigma} \big( \mathds{P}_\tau^{N_\tau}[{\bm t}_1] -  \mathds{P}_\tau^{N_\tau}[{\bm t}_2] \big) + \big( {\bm \gamma}^\phi \big)^T \mathbb{h_\sigma} \big(  \mathds{P}_\tau^{N_\tau}[{\bm \phi}_1] -  \mathds{P}_\tau^{N_\tau}[{\bm \phi}_2] \big)\\
\nonumber+& \big( \tilde{\bm \gamma}^t \big)^T \mathbb{h_\sigma} \big(  \mathds{P}_\tau^{N_\tau}[(\mathds{D}_\tau{\bm t}_1)] -  \mathds{P}_\tau^{N_\tau}[(\mathds{D}_\tau{\bm t}_2)] \big)\\
\nonumber+& \big( \tilde{\bm \gamma}^\phi \big)^T \mathbb{h_\sigma} \big(  \mathds{P}_\tau^{N_\tau}[(\mathds{D}_\tau{\bm \phi}_1)] -  \mathds{P}_\tau^{N_\tau}[(\mathds{D}_\tau{\bm \phi}_2)] \big)\\
\nonumber+& \big( {\bm \kappa}^\phi \big)^T \mathbb{h_\tau} \big(  \mathds{P}_\sigma^{0}[{\bm \phi}_1] - {\bm 0} \big) + \big( \tilde {\bm \kappa}^\phi \big)^T \mathbb{h_\tau} \big(  \mathds{P}_\sigma^{N_\sigma}[{\bm \phi}_1] - {\bm 0} \big)\\
+& \big( {\bm \xi}^\phi \big)^T \mathbb{h_\tau} \big(  \mathds{P}_\sigma^{0}[{\bm \phi}_2] - {\bm 0} \big) + \big( \tilde {\bm \xi}^\phi \big)^T \mathbb{h_\tau} \big(  \mathds{P}_\sigma^{N_\sigma}[{\bm \phi}_2] - {\bm 0} \big).\label{eq:discrEsim}
\end{align}

For the discrete Noether charge corresponding to <ref>, we must implement spatial integration as described in <ref>. To this end we construct the $N_\tau\times (N_\tau \cdot N_\sigma)$ matrix $\mathds{H}_\sigma$. Applying $\mathds{H}_\sigma$ to a discretized function $\bm f$ results in an array of $N_\tau$ length. Here the $N_\sigma$ component vector $({\bm h}_\sigma)_i=(\mathbb{h}_\sigma)_{ii}$, which contains the diagonal entries of the spatial quadrature matrix $\mathbb{h}_\sigma$. The matrix $\mathds{H}_\sigma$ then has the form (<ref>).

To obtain the array for the Noether charge with $N_\tau$ entries we thus apply

\begin{align}
\nonumber {\bm Q}_t^{\rm L}= \mathds{H}_\sigma &\underbracket{\Big\{ ( \mathds{D}_\tau{\bm t_1}) + \frac{1}{T} \Big( (\mathds{D}_\sigma {\bm \phi}_1)^2\circ(\mathds{D}_\tau{\bm t_1)}- (\mathds{D}_\tau {\bm \phi}_1)\circ(\mathds{D}_\sigma {\bm \phi}_1)\circ(\mathds{D}_\sigma {\bm t}_1)\Big)\Big\}}_{{\bm J}^0\in {\mathbb{R}^{N_\tau\times N_\sigma}}}\\
+  & \underbracket{\Big\{ ({\bm h}_\sigma^T \tilde{\bm \lambda}^t)  \mathfrak{d}^\tau[0]+ ( {\bm h}_\sigma^T \tilde{\bm \gamma}^t ) \mathfrak{d}^\tau[N_\tau]\Big\}}_{\rm Lagr.\, mult.\, contrib.},\label{eq:dsicrNoetherl}
\end{align}

where $\mathfrak{d}^\tau[k]$ refers to an array with $N_\tau$ entries, defined in (<ref>).

To confirm quantitatively the conservation of this Noether charge, let us introduce $\Delta {\bm E}= {\bm Q}_t^{\rm L} - {\bm Q}_t^{\rm L}[0]$ the difference between the discrete Noether charge ${\bm Q}_t^{\rm L}$ and its initial value, which, due to its dependence on only the derivatives of the degrees of freedom, is fully specified by the initial conditions.

As we saw for the Noether charge, the presence of Lagrange multiplier terms, implementing the initial and connection conditions, affects the values of our system at the final (and potentially initial) time slice.
The influence of the Lagrange multipliers manifests itself in the solutions of both the coordinates and fields, which show deviations from the discretized equations of motion on the boundary (for a detailed discussion see ref. [44]).
Note that in the interior the solutions we obtain exactly fulfill the discretized equations of motion and the only deviations occur on the boundary. We will explicitly demonstrate that these deviations on the boundary do not affect the global convergence to the correct solution under grid refinement.

§.§ Numerical results

We next present numerical results that exemplify how our novel approach addresses two of the three central challenges to formulating and discretizing IBVPs

* Due to the presence of a manifest continuum symmetry (time translations) after discretization, the discrete Noether charge ${\bf Q}_t^{\rm L}$ remains preserved over the whole simulation period.

Since our SBP operators exactly mimic integration by parts, the value of the Noether charge remains conserved at its initial value set solely by the initial data.

* The dynamical coordinate map $t(\tau,\sigma)$ is not a trivial sheet but exhibits different values of $\tau$ and $\sigma$ derivatives throughout the simulation domain. This amounts to regions of finer and coarser temporal resolution, constituting a dynamically emerging adaptive mesh, i.e. an automatic AMR procedure.

Our simulation domain is defined by the parameters $\tau\in[0,1/2]$ and $\sigma\in[0,1]$.
Choosing $N_\tau=60$ and $N_\sigma=48$ for the presentation of the main results below, we have $\Delta \tau=1/118$ and $\Delta \sigma =1/47$. Our choice of grid points and spacings represents a compromise between capturing relevant wave dynamics and computational cost.

The physical extent of our simulation domain is determined by the choice of coordinate mappings. As we use the trivial mapping $x=\sigma$ in spatial direction, the total x-direction extent is equal to that of $\sigma$. On the other hand the temporal mapping is fully dynamic and we control its extent by prescribing an initial value $t(\tau^{\rm i},\sigma)=0$ and velocity of $\td(\tau^{\rm i},\sigma)=5/2$.

We choose a value of $T=10^4$ such that deviations from the simple sheet behavior of the time mapping will be of at most on the percent level. In realistic systems, one will encounter an even larger value of $T$, since the corrections to the behavior encoded in the conventional field theory action $S_{\rm ft}$ in <ref> are minute. What is important is that $T$ is present, so that the coordinate maps can be considered as dynamical. From a practical point of view, as long as $T$ is chosen much larger than the energy density of the field, encoded in the initial conditions, the solution obtained for the field will benefit from the automatic AMR procedure and also converge to a continuum result that is indistinguishable from that of the conventional action.

The dynamical adaptation of the time-mapping leads the simulation to end at a value of the time coordinate close to $t^{\rm f}\approx5/4$. We choose this value, as our wave packages propagate with physical velocity $v=dx/dt=1$. Between $t^{\rm i}$ and $t^{\rm f}$ the wave will have reflected from the hard spatial system boundaries, characterized by $\phi=0$, and will have met again at the center of the simulation domain. This non-trivial dynamics allows us to showcase relevant aspects of our simulation approach.

We initialize the field $\phi$ with a localized bump, which is symmetrically located in the center of the spatial domain at $\sigma=1/2$ with functional form $\phi(0,\sigma)={\rm sin}[\pi\sigma]{\rm exp}[-100(\sigma-1/2)^2]$. On the spatial boundary we prescribe vanishing Dirichlet boundary conditions. In <ref> we plot the initial and boundary values for the dynamical field $\phi$ as red points and the initial value of $t$ as blue points.

Initial and boundary values provided for the dynamic field $\phi$ as red points and initial values of the time mapping $t$.

In order to mimic the continuum theory, we deploy SBP operators. The examples shown will be based on the  approximation of the derivative. At the end of the section we will provide some comparisons with results based on  operators.

(left) Eigenvalue spectrum of the unregularized  finite difference operators $\mathds{D}_\tau$ (red circles) and $\mathds{D}_\sigma$ (blue crosses) on a grid with $N_\sigma=24$ and $N_\tau=16$. Since $\Delta \tau<\Delta \sigma$ the purely imaginary eigenvalues of $\mathds{D}_\tau$ spread over a larger interval than for the spatial derivative. Each operator has exactly two zero modes depicted at the origin. (right) Eigenvalue spectrum of the regularized  operators (using the initial and boundary data from <ref>) in affine coordinates. The zero modes of the SBP operator are lifted and the original physical zero mode, i.e. the constant function, is now associated with the unit eigenvalue.

An important ingredient for use of SBP operators in the discretization of IBVP actions is their regularization, in order to avoid the $\pi$-mode. In  <ref> we provide for completeness the actual eigenvalue spectrum of our  operators for a smaller grid with $N_\sigma=24$ and $N_\tau=16$ points, which already allows us to showcase all relevant properties. The left panel shows the purely imaginary eigenvalues of $\mathds{D}_\tau$ (red circles) and $\mathds{D}_\sigma$ (blue crosses) with their two degenerate zero eigenvalues located at the origin. In the right panel, we instead show the spectrum of the three regularized operators deployed in <ref> each formulated in affine coordinates with a structure as shown in <ref>. We have two temporal ones $\bar{\mathds{D}}^t_\tau$ and $\bar{\mathds{D}}^\phi_\tau$, as well as the spatial one $\bar{\mathds{D}}^\phi_\sigma$. The data used for the regularization are the initial and spatial boundary values shown in <ref>. Note that all zero modes have been successfully lifted and that the original physical zero mode, the constant function, is now associated with an eigenvalue of unit magnitude.

In order to determine the classical solution we carry out a numerical optimization procedure, locating the critical point of <ref>
in the full set of dynamical degrees of freedom ${\bf d}=\{ {\bf t}_{1,2},{\bm \phi}_{1,2}, {\bm \lambda}^t,{\bm \lambda}^\phi,\tilde {\bm \lambda}^t,\tilde{\bm \lambda}^\phi,{\bm \gamma}^t,\tilde {\bm \gamma}^\phi,\tilde {\bm \gamma}^t,{\bm \gamma}^\phi, {\bm \kappa}^\phi, \tilde{\kappa}^\phi,{\bm \xi}^\phi,\tilde{\bm \xi}^\phi\}$.

For our proof-of-principle study we use the   software[Example  code, based on the  operator is available at the Zenodo repository under open access [52].] and the minimization routines it provides. Note that in the presence of Lagrange multipliers the critical point of the action functional is in general a saddle point and not an extremum. With most established optimization algorithms designed to locate extrema, we first convert the action functional into a suitable form, i.e. we optimize not on $\mathds{E}^{\rm L}_{\rm IBVP}$ but on the norm of its gradient $|{\bm \nabla}_{\bm d}\mathds{E}^{\rm L}_{\rm IBVP}|$. When $\mathds{E}^{\rm L}_{\rm IBVP}$ features a saddle point, then the norm of its gradient exhibits a genuine minimum, since its values must be positive definite. In this way we convert saddle points into extrema, which can be straight forwardly located.

In practice we use multi-step preconditioning to locate the solution. First, in order to obtain appropriate starting values, we optimize on an action where the coordinate maps and field degrees are decoupled, similar to the $T\to\infty$ limit of the full theory. To exploit optimally the different stability and convergence properties of optimization algorithms, we start with the  method of the  command , iterating for around one to two-thousand steps, followed by application of the interior-point optimization [53] based method . With this solution in hand we turn to the full action and use the  method to obtain the final result.

Since the  method in  makes reference to an external compiled library, it is unable to benefit from the inbuilt high precision arithmetic, which limits us to find extrema up to residuals of $\Delta\sim 10^{-20}$. While at first this may appear more than satisfactory, it was observed in the study of the world-line formalism in [38] that even smaller global residuals are often needed to e.g. enforce the physical limit, i.e. the equality of the degrees of freedom on the forward and backward branch, down to zero in double precision. The reason for the occurrence of such small residuals lies in the intricate cancellations that take place between the degrees of freedom on the forward and backward branch. Hence, as we are limited to using , our overall precision will be around $\Delta\sim10^{-6}$. I.e. within the accuracy of our minimization
a numerical zero amounts to values smaller than $\Delta\sim10^{-6}$.

We proceed by presenting the first result of our numerical study. The top panel of <ref> depicts the solution obtained for the field $\phi(\tau,\sigma)$ as a function of our abstract parameters $\tau$ and $\sigma$. The bottom panel on the other hand contains the values for the time-mapping $t(\tau,\sigma)$, which, due to our choice of $T=10^4$ at first sight resembles a plain sheet.

(top) Two  Classical solution for the field $\phi$ (top) and the time-mapping $t$ (bottom) from the critical point of <ref> evaluated on a grid with $N_\sigma=48$ and $N_\tau=60$ points. Note that the initial and boundary conditions of <ref> are exactly obeyed.

Our initial field configuration amounts to a  spatially symmetric bump of nearly Gaussian form. Based on intuition, we expect that this structure corresponds to an almost equal amount of left and right propagating modes. Indeed, we see two wave-packages emerge from the initial bump, travelling with opposite but equal speed towards the spatial boundaries of the simulation domain.

The vanishing Dirichlet boundary conditions amount to the physical situation of the wave medium (such as e.g. an elastic sheet) being fastened at the boundary. As one would expect from Newton's third law, the amplitude of the wave-package changes sign after reflection. The wave attempts to push the boundary upward and the
boundary exerts an equal but opposite force onto the medium, leading instead to a downward movement. After reflection the wave-packages travel back towards the center of the spatial domain, where they interfere before heading towards the opposite end of the spatial domain.

Next, we take a closer look at the behavior of the time mapping, whose $\tau$ derivative is plotted in the top panel of <ref>, while its $\sigma$ derivative is shown in the bottom panel. We observe nontrivial structures intimately related to the dynamics of the propagating field.

The temporal (top) and the spatial (bottom) derivatives of the time-mapping $t(\tau,\sigma)$ as obtained from the critical point of <ref>, evaluated on a grid with $N_\sigma=48$ and $N_\tau=60$ points.

Starting with the $\tau$ derivative, if coordinate maps and fields evolved independently (as would be the case for $T\to\infty$), one would expect to find a constant value of $\Delta t/\Delta \tau=\dot t_{\rm IC}=5/2$ . This constant value amounts to a uniform temporal grid spacing throughout the simulation domain. For finite values of $T$ we observe deviations from such constant behavior, with larger values indicating a coarser and smaller values indicating a finer temporal grid in different regions of the domain spanned by $(\tau,\sigma)$. In particular two structures are noteworthy: close to the center of the $\sigma$ domain we find a double ridge structure stretching out over a finite interval in $\tau$ direction and then repeating as $\tau$ progresses. In addition we find significant deviations from constancy at the spatial boundaries of the simulation domain.