Holger Heumann
email: holger.heumann@inria.fr

Ralf Hiptmair
email: hiptmair@sam.math.ethz.ch

Cecilia Pagliantini
email: cecilia.pagliantini@sam.math.ethz.ch

STABILIZED GALERKIN FOR TRANSIENT ADVECTION OF DIFFERENTIAL FORMS

Keywords: 2010 Mathematics Subject Classification. Primary: 65M60, 65M12; Secondary: 65L06 Lie derivative, Discrete di↵erential forms, Stabilized Galerkin, explicit Runge Kutta timestepping

We deal with the discretization of generalized transient advection problems for di↵erential forms on bounded spatial domains. We pursue an Eulerian method of lines approach with explicit timestepping. Concerning spatial discretization we extend the jump stabilized Galerkin discretization proposed in [H. Heumann and R. Hiptmair, Stabilized Galerkin methods for magnetic advection, Math. Modelling Numer. Analysis, 47 (2013Analysis, 47 ( ), pp. 1713Analysis, 47 ( -1732] ] to forms of any degree and, in particular, advection velocities that may have discontinuities resolved by the mesh. A rigorous a priori convergence theory is established for Lipschitz continuous velocities, conforming meshes and standard finite element spaces of discrete di↵erential forms. However, numerical experiments furnish evidence of the good performance of the new method also in the presence of jumps of the advection velocity.

1. Introduction. The equations of magneto-hydrodynamics (MHD) [START_REF] Davidson | An Introduction to Magnetohydrodynamics[END_REF]Section 3.8], [START_REF] Goedbloed | Principles of Magnetohydrodynamics[END_REF]Section 4.1] provide a consistent description of the interaction of electromagnetic fields and conducting non-magnetic fluids like plasmas. The standard model for resistive MHD under a quasi-neutrality assumption comprises balance equations for mass, momentum and energy together with material laws and Maxwell's equations in their magneto-quasistatic reduction (eddy current model) for the electromagnetic fields.

The traditional formulation of the linear eddy current model in the presence of a conducting fluid moving with velocity = (x, t) boils down to the evolution PDE

@ t u + curl"curl u + ↵u + curl u ⇥ + grad( • u) = f , (1)
governing the evolution of the unknown magnetic vector potential u and with " being the magnetic di↵usion coe cient. Alternatively, one may rely on the magnetic induction field as unknown, again denoted by u, which yields @ t u grad"divu + ↵u + divu + curl(u ⇥ ) = f .

(2)

1.1. Generalized advection-di↵usion evolution problem. Both (1) and [START_REF] Ambrosio | The flow associated to weakly di↵erentiable vector fields: recent results and open problems[END_REF] belong to a single family of second order evolution problems, which we have dubbed generalized advection-di↵usion problems. For a unified statement we rely on the language of exterior calculus. In that notation, the strong form of the generalized advection-di↵usion equation in the space-time domain ⌦ ⇥ I, I := [0, T ], written in terms of di↵erential forms reads ?@ t !(t) + ( 1) k+1 d" ? d!(t) + ?↵!(t) + ?L !(t) = f (t), in ⌦ ⇥ I, tr !(t) = tr g(t), on ( in [ 0 ) ⇥ I, tr(i n !(t)) = tr s(t), on in ⇥ I,

!(0) = ! 0 , in ⌦. (3)
where !(t) is a time-dependent di↵erential k-form on the bounded domain ⌦ ⇢ R n , : ⌦ ⇥ I ! R n is a given velocity field and f (t) 2 ⇤ n k (⌦) a source term. The scalar di↵usivity parameter " and the reaction coe cient ↵ are non-negative and bounded functions ⌦ ! R, and the boundary conditions are imposed at the inflow boundary in := {x 2 @⌦ : • n(x) < 0} and at the "elliptic boundary" 0 (where the di↵usion parameter " > 0). All other notations are borrowed from [25, Section 2] and summarized in Table 1.

Symbol

Meaning in exterior calculus ⇤ k (⌦) : space of (smooth) di↵erential forms on a bounded domain ⌦ ⇢ R n d : exterior derivative operator ⇤ k (⌦) ! ⇤ k+1 (⌦) [ The so-called vector proxies1 establish the connection between (1), ( 2) and (3). Indeed, in R 3 endowed with the Euclidean inner product natural isomorphisms between ⇤ k (R 3 ) and R or R 3 can be defined. The fields associated to di↵erential forms are called proxy fields for the forms and exterior calculus operations on forms correspond to operations on scalar functions and vector fields, see [START_REF] Heumann | Convergence of Lowest Order Semi-Lagrangian Schemes[END_REF]Section 2], [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF]Table 2.1], or Table 2 andTable 3. From these identifications we see that (1), (2) correspond to (3) for k = 1 and k = 2, respectively.

! 2 ⇤ k (⌦) k = 0 k = 1 k = 2 k = 3 d! gradu curl u divu i ! • u u ⇥ u ! divu curl u gradu j ! u u ⇥ • u L ! • gradu grad( • u)+curl u⇥ curl (u ⇥ ) + divu div(u ) L ! div(u ) curl ( ⇥ u) divu ⇥curl u grad( • u) • gradu tr u(x) n(x) ⇥ u(x) u(x) • n(x) H⇤ k (⌦) H 1 (⌦) H(curl, ⌦) H(div, ⌦) L 2 (⌦)
^: ⇤ 1 (R 3 ) ⇥ ⇤ 1 (R 3 ) ! ⇤ 2 (R 3 ) ⇥ : R 3 ⇥ R 3 ! R 3 (cross product) ^: ⇤ 1 (R 3 ) ⇥ ⇤ 2 (R 3 ) ! ⇤ 3 (R 3 ) • : R 3 ⇥ R 3 ! R (dot product)
Hodge operator ? ?

: ⇤ 0 (R 3 ) ! ⇤ 3 (R 3 ) id: R ! R ? : ⇤ 1 (R 3 ) ! ⇤ 2 (R 3 ) id: R 3 ! R 3
Table 3. Correspondence between wedge product and Hodge operator for di↵erential forms and vector proxies. More details in [27, Section 2.2], [5, Table 2.1].

For k = 0 the evolution operator in (3) written in vector proxies becomes the familiar and widely studied second order advection-di↵usion equation for the unknown scalar function u @ t u div"gradu + ↵u + • gradu = f.

(4)

By analogy we conclude that in the generalized advection-di↵usion problem (3), the di↵usion operator is d ? d, the zero-th order term amounts to a reaction term and the advection operator is the Lie derivative L associated to the velocity field .

It is well known that for scalar advection-di↵usion equation ( 4) straightforward Galerkin discretization with Lagrangian finite elements will break down in the singular perturbation limit of vanishing di↵usion. Thus, robustness for " & 0 will also be a key issue for the spatial discretization of (1) and [START_REF] Ambrosio | The flow associated to weakly di↵erentiable vector fields: recent results and open problems[END_REF]. In this article we tackle the challenge of robust Eulerian spatial finite element discretization for the general advection-di↵usion problem [START_REF] Arnold | Spaces of finite element di↵erential forms[END_REF]. In fact, we will focus on the pure advection problem obtained from (3) for " = 0; if a scheme performs well in this case, it will also be suitable for (3) when augmented with a standard Galerkin discretization of the di↵usion term.

1.2. Pure advection problem: Statement and well-posedness. We introduce the spaces

V := {! 2 L 2 ⇤ k (⌦) : L ! 2 L 2 ⇤ k (⌦) , Z in tr i (! ^?!) < 1}, W := {! 2 V : tr ! = g, tr i n ! = s on in , g(t) 2 L 2 ⇤ k ( in ), s(t) 2 L 2 ⇤ k 1 ( in )},
and state the pure advection initial boundary value problem: For ! 0 2 W |t=0 and

f 2 C 0 (I; L 2 ⇤ k (⌦)) find ! 2 C 1 (I; L 2 ⇤ k (⌦)) \ C 0 (I; W ) such that @ t !(t) + ↵!(t) + L !(t) = f (t), in ⌦ ⇥ I , !(0) = ! 0 , in ⌦ . (5)
If does not depend on t and 2 W 1,1 (⌦), the Hille-Yosida theorem in [19, Theorem 6.52] can be directly applied to show that the variational problem associated to (5

): find ! 2 C 1 (I; L 2 ⇤ k (⌦)) \ C 0 (I; W ) such that, for all ⌘ 2 L 2 ⇤ k (⌦) (@ t !, ⌘) ⌦ + (↵!, ⌘) ⌦ + (L !, ⌘) ⌦ = (f, ⌘) ⌦ (!(0), ⌘) ⌦ = (! 0 , ⌘) ⌦ is well-posed. Here (•, •) ⌦ denotes the L 2 ⇤ k (⌦) inner product (!, ⌘) ⌦ := R ⌦ ! ^?⌘.
Further, for velocity fields uniformly continuous in time and Lipschitz continuous in space, that is, 2 C 0 (I; W 1,1 (⌦)), it can be shown [START_REF] Heumann | Convergence of Lowest Order Semi-Lagrangian Schemes[END_REF]Lemma 3.4] that the monotonicity condition

Z ⌦ ↵ + 1 2 (L (•,t) + L (•,t) ) ! ^?! ↵ 0 Z ⌦ ! ^?! 8! 2 L 2 ⇤ k (⌦) , 8t 2 I , (6)
for some constant ↵ 0 > 0 and with L = ( 1) k(n k) ? L ?, ensures that the operator ↵ id +L : [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] falls into the class of evolution problems for the so-called Friedrichs' symmetric operators [START_REF] Friedrichs | Symmetric positive linear di↵erential equations[END_REF] and then [30, pp. 143-145] gives well-posedness of (5

W ! L 2 ⇤ k (⌦)
) if ⌦ = R n .
These results require to be Lipschitz continuous in space. However, MHD solutions feature shocks that give rise to discontinuous velocities; discontinuous transport velocities are relevant in the context of magneto-quasistatic Maxwell's equations, also in the limit of small di↵usion.

A well-posedness theory for velocity fields with less regularity is available only for scalar advection. In [START_REF] Diperna | Ordinary di↵erential equations, transport theory and Sobolev spaces[END_REF] DiPerna and Lions showed well-posedness of the scalar advection problem for velocity fields 2 L 1 loc (0, T ;

W 1,1 (R n )) with div 2 L 1 (0, T ; L 1 (R n ))
through the concept of renormalized solutions. More recently, Ambrosio in [START_REF] Ambrosio | Transport equation and Cauchy problem for BV vector fields[END_REF] provided an extension of this breakthrough to transport velocity fields in L 1 loc (0, T ; BV loc (R n )) and div 2 L 1 (0, T ; L 1 loc (R n )). Moreover, a notion of generalized flow associated with low regular velocity fields (the regular Lagrangian flow) and an extension of the characteristics theory to beyond the smooth context have been subject of investigation of several authors, see [START_REF] Crippa | Regularity and compactness for the DiPerna-Lions flow[END_REF], [START_REF] Ambrosio | The flow associated to weakly di↵erentiable vector fields: recent results and open problems[END_REF], [START_REF] Bouchut | Lagrangian flows for vector fields with gradient given by a singular integral[END_REF] and the references therein. To the best of our knowledge, beside the case of scalar transport, a well-posedness theory for the generalized transport problem [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] with low regular advection velocities has not been developed.

Even though the above mentioned results have been established for nearly incompressible velocity fields (see [START_REF] Lellis | Notes on hyperbolic systems of conservation laws and transport equations[END_REF] for a detailed overview), the assumption on the boundedness of the divergence of the velocity (absolute continuity with respect to the Lebesgue measure in the BV case) is of crucial importance for the well-posedness of the scalar advection problem. In the context of the generalized transport problem for a di↵erential k-form, this corresponds to require the operator L + L to be bounded in space, which conceals a rather strong assumption on the regularity of the velocity itself, when k = 1, 2.

1.3. Novelty and outline. A full discretization of (5) was already presented in [START_REF] Heumann | Convergence of Lowest Order Semi-Lagrangian Schemes[END_REF]. There, the authors introduced a semi-Lagrangian approach. Conversely, in the present paper we pursue a mesh based Eulerian method of lines approach to (5), employing a (jump) stabilized Galerkin discretization and piecewise polynomial discrete di↵erential forms for spatial discretization. Our new methods will be constructed to accommodate discontinuous velocities aligned with the mesh.

A jump-stabilized discontinuous Galerkin method for the stationary advection problem for 0-forms in R 3 and Lipschitz continuous velocities 2 W 1,1 (⌦), was introduced and theoretically analyzed in [START_REF] Brezzi | Discontinuous Galerkin methods for first-order hyperbolic problems[END_REF]. An extension of these results to the magnetic advection problem (1-forms in R 3 , cf. (1)) was proposed in [START_REF] Heumann | Stabilized Galerkin methods for magnetic advection[END_REF], where a priori convergence rates were derived for both fully discontinuous piecewise polynomial functions and H(curl, ⌦)-conforming finite elements. Discontinuous velocity fields were not taken into account. We remark that for discontinuous velocities, even the spatial discretization of the scalar transport problem (4), for which existence and uniqueness of weak solutions are known, is discussed only rarely ( [START_REF] Boyer | Analysis of the upwind finite volume method for general initial-and boundary-value transport problems[END_REF], [START_REF] Walkington | Convergence of the discontinuous Galerkin method for discontinuous solutions[END_REF]).

The remainder of the paper is organized as follows. In Section 2, we devise a stabilized Galerkin spatial semi-discretization for the generalized stationary advection problem (5) for merely piecewise smooth velocity . It is an extension of the method introduced in [START_REF] Brezzi | Discontinuous Galerkin methods for first-order hyperbolic problems[END_REF] for 0-forms and in [START_REF] Heumann | Stabilized Galerkin methods for magnetic advection[END_REF] for 1-forms. Trial and test functions are polynomial discrete di↵erential forms, which will be introduced in Section 2.5. The new method is a substantial extension of the scheme presented in [START_REF] Heumann | Stabilized Galerkin methods for magnetic advection[END_REF] to forms of arbitrary degree, any spatial dimension and velocities with jumps.

Next, Section 3 establishes stability a priori convergence estimates for the stabilized Galerkin discretization in the stationary setting. For want of well-posedness results for the generalized advection problem in case of discontinuous , these investigations are confined to Lipschitz continuous velocities 2 W 1,1 (⌦). The stability and consistency results obtained in that section are instrumental for the convergence analysis of the fully discrete scheme in Section 4. We study explicit timestepping following the approach of [START_REF] Levy | From semidiscrete to fully discrete: stability of Runge-Kutta schemes by the energy method[END_REF] and [START_REF] Burman | Explicit Runge-Kutta schemes and finite elements with symmetric stabilization for first-order linear PDE systems[END_REF].

Finally, in Section 5 and Section 6 the performance of the new method is tested in various numerical experiments for both the stationary and transient generalized advection problem (5) in 2D. The tests cover both continuous and discontinuous velocities and employ tensor product grids and triangular meshes.

2. Spatial discretization.

2.1.

Stationary generalized advection problem. The Eulerian method of lines policy applies timestepping after discretization in space. Therefore, we will first address the spatial discretization of (5) and we start from the stationary generalized advection boundary value problem for a k-form ! on the bounded computational domain ⌦ ⇢ R n :

! + L ! = f, in ⌦ , (7a)
tr ! = g, on in , (7b)

tr i n ! = s, on in , (7c)
with f 2 L 2 ⇤ k (⌦), g 2 L 2 ⇤ k ( in ), s 2 L 2 ⇤ k 1 ( in )
, and piecewise Lipschitz continuous velocity field . As stated in [23, p. 59], if 2 W 1,1 (⌦) problem ( 7) is well-posed in V under the assumption (6).

2.2. Transmission conditions. We aim for stabilized Galerkin methods that, crudely speaking, involve a penalization of suitable jumps across interfaces inside ⌦. In order to select the right jump terms, we have to understand the natural transmission conditions across an internal interface f ⇢ ⌦ satisfied by a solution ! of [START_REF] Bossavit | On the geometry of electromagnetism. (2): 'Geometrical objects[END_REF].

For smooth velocity

2 W 1,1 (⌦), the requirement L ! 2 L 2 ⇤ k (⌦) read in distributional sense, involves the transmission condition tr [i (! ^?⌘)] f = 0 8 ⌘ 2 C 1 0 ⇤ k (⌦), (8)
for any oriented (piecewise) smooth n 1-dimensional surface f ⇢ ⌦. Here, we wrote [•] f for the jump of a function across the surface f . This formula is a consequence of the integration by parts formula for the Lie derivative

(L !, ⌘) ⌦ (!, L ⌘) ⌦ = Z @⌦ tr i (! ^?⌘) 8!, ⌘ 2 C 1 ⇤ k (⌦). (9)
The transmission conditions (8) carry over to Lipschitz continuous velocity 2 W 1,1 (⌦). Clearly, no transmission conditions are imposed across surfaces tangential to (characteristic surfaces).

In case of discontinuous velocity , an interpretation of L ! in the sense of distributions is no longer available. Therefore, at jumps of resort to a strong interpretation of L !. Appealing to Cartan's homotopy formula (see for example [START_REF] Schwarz | Hodge Decomposition-a Method for Solving Boundary Value Problems[END_REF]Equation 2.3] or [START_REF] Lee | Introduction to smooth manifolds[END_REF]Theorem 14.35]) L = di + i d, we conclude the strong transmission conditions

tr [!] f = tr [i !] f = 0 8 oriented surfaces f ⇢ ⌦, [ ] f 6 = 0, ( 10
) from demanding ! 2 L 2 ⇤ k (⌦), i ! 2 L 2 ⇤ k 1 (⌦) and d! 2 L 2 ⇤ k+1 (⌦).
2.3. Stabilized Galerkin variational formulation. In the following, let T h = {T } be a cellular partition (generalized triangulation) of ⌦ ⇢ R n into (curved) polyhedra T . Denote by F and F @ the set of interior and boundary n 1-faces of T h (named facets) and F = F [ F @ . The set of facets at the inflow boundary is defined as F @ := {f 2 F @ : f ⇢ in } and in = [ f 2F @ f , whereas F @ \ F @ is the set of facets at the outflow boundary. An oriented facet f has a distinguished normal n f . Any facet f , as part of the boundary of some element T 2 T h , has either

n f = n T | f or n f = n T | f . Then, given ! 2 ⇤ k (⌦)
, its two di↵erent restrictions to f are denoted by ! + and ! , e.g. ! + := ! | T + where element T + has outward normal n f . Hence, we can introduce the notion of jump and average across a facet

f 2 F as [!] f := ! + ! , {!} f := 1 2 (! + + ! ).
For f ⇢ @⌦ we assume f to be oriented such that n f points outwards and [!] f = {!} f := !. We also write h T := diam T and h := max T 2T h h T .

Further, let ⇤ k h (T h ) denote some piecewise polynomial approximation space for di↵erential k-forms. Here ⇤ k h (T h ) could be either a H⇤ k (⌦)-conforming space

⇤ k h (T h ) ⇢ H⇤ k (⌦) or a non-conforming space ⇤ k h (T h ) ⇢ L 2 ⇤ k (⌦) for which ⇤ k h (T h ) 6 ⇢ H⇤ k (⌦) .
The method is formulated in the general framework of time-dependent velocity fields = (x, t) and relies on the assumption that the possible (space) discontinuities of the velocity are resolved by the mesh: Assumption 1. For every t 2 I we have (•, t) | T 2 W 1,1 (T ) for each T 2 T h , that is the velocity field is assumed to be T h -piecewise Lipschitz continuous.

This may seem to be a severe limitation but for our purposes it represents a reasonable condition in view of the fact that the velocity field is obtained from numerically solving the MHD system.

Next, multiplying equation (7a) by a test form ⌘ h 2 ⇤ k h (T h ) and applying the integration by parts rule [START_REF] Boyer | Analysis of the upwind finite volume method for general initial-and boundary-value transport problems[END_REF], results in

(↵! h , ⌘ h ) ⌦ + X T 2T h (! h , L ⌘ h ) T + X T 2T h Z @T tr i (! h ^?⌘ h ) = (f, ⌘ h ) ⌦ 8⌘ h 2 ⇤ k h (T h ) .
Let j be the formal adjoint of the contraction operator i as in Table 1. Applying the following product rule

i (! ^?⌘) = i ! ^?⌘ + ( 1) k+`! ^?j ⌘ 8 ! 2 ⇤ k (⌦) , ⌘ 2 ⇤ `(⌦) (11)
to the boundary terms, results in

(↵! h , ⌘ h ) ⌦ + X T 2T h (! h , L ⌘ h ) T + X f 2F Z f tr [i ! h ^?⌘ h ] f + X f 2F Z f tr [! h ^?j ⌘ h ] f = (f, ⌘ h ) ⌦ 8 ⌘ h 2 ⇤ k h (T h ) .
Moreover, it can be easily verified that, for all

µ h , ⌘ h 2 ⇤ k h (T h ), it holds X f 2F Z f tr [µ h ^?⌘ h ] f = X f 2F Z f tr({µ h } f ^?[⌘ h ] f ) + X f 2F Z f tr([µ h ] f ^?{⌘ h } f ).
For ! 2 W solution of problem [START_REF] Bossavit | On the geometry of electromagnetism. (2): 'Geometrical objects[END_REF], the transmission conditions (10) at the mesh facets

tr [!] f = tr [i !] f = 0 8 f 2 F , yield the variational formulation: find ! h 2 ⇤ k h (T h ) such that a h (! h , ⌘ h ) = l(⌘ h ) for all ⌘ h 2 ⇤ k h (T h ), where l(⌘ h ) := (f, ⌘ h ) ⌦ X f 2F @ Z f tr i (g ^?⌘ h ) ( 12
)
a h (! h , ⌘ h ) := (↵! h , ⌘ h ) ⌦ + X T 2T h (! h , L ⌘ h ) T + X f 2F @ \F @ Z f tr i (! h ^?⌘ h ) + X f 2F Z f tr({i ! h } f ^?[⌘ h ] f ) + Z f tr({! h } f ^?[j ⌘ h ] f ). ( 13
)
As it is well-known, classical Galerkin finite element discretization of advection problems su↵er from instabilities. Therefore, devising stabilization techniques to counteract this limitation has been investigated widely. We consider the following stabilization operator, for all

⌘ h 2 ⇤ k h (T h ), s h (! h , ⌘ h ) := X f 2F Z f c f tr([i ! h ] f ^?[⌘ h ] f ) + Z f cf tr([! h ] f ^?[j ⌘ h ] f ), (14)
where the stabilization scalar functions c f (x) and cf (x) may depend on the velocity field and on the facets diameter h f . Throughout, the stabilization parameters are assumed to satisfy the following:

Assumption 2. We assume that c f (x) and cf (x) satisfy:

c f • n f c 0 > 0 and cf • n f c0 > 0 uniformly for all facets f 2 F .
In particular, by considering the direction of the numerical fluxes as given by the average of the velocity field, the choice

c f = cf = 1 2 { } f • n f |{ } f • n f | , f 2 F , (15)
gives a scheme with upwind fluxes (see [START_REF] Heumann | Eulerian and Semi-Lagrangian Methods for Advection-Di↵usion of Di↵erential Forms[END_REF]Remark 4.1.2] in the case 2 W 1,1 (⌦)). Indeed, from (13) together with ( 14) the facets contribution, for

! h , ⌘ h 2 ⇤ k h (T h ), reads X f 2F Z f tr({i ! h } f ^?[⌘ h ] f ) + Z f c f tr([i ! h ] f ^?[⌘ h ] f ) + X f 2F Z f tr({! h } f ^?[j ⌘ h ] f ) + Z f c f tr([! h ] f ^?[j ⌘ h ] f ) = 1 2 X f 2F Z f tr ⇣ (1 + 2c f )(i ! h ) + ^?[⌘ h ] f ⌘ + tr ⇣ (1 2c f )(i ! h ) ^?[⌘ h ] f ⌘ + 1 2 X f 2F Z f tr ⇣ (1 + 2c f )! + h ^?[j ⌘ h ] f ⌘ + tr ⇣ (1 2c f )! h ^?[j ⌘ h ] f ⌘ .
Note that, since the velocity field is discontinuous, the upwind direction at the mesh facets may not be well defined. Here we consider the direction of the stream as the one given by the average of the velocity. However, other possibilities are feasible: an upwind direction given locally by the velocity field can be used, even if this choice will lead to non-unique numerical fluxes at mesh facets.

The evaluation of the terms in (13) involving the Lie derivative L ⌘ h requires the knowledge of the first order derivatives of the velocity field . Note that since the velocity is assumed to be a smooth function in all elements T 2 T h , the quantity (! h , L ⌘ h ) T is well defined for all T 2 T h . However, as suggested in [START_REF] Heumann | Stabilized Galerkin methods for magnetic advection[END_REF], a di↵erent equivalent formulation of the bilinear form a h (•, •) is convenient for implementation purposes.

Proposition 1. The following equality holds for all

! h , ⌘ h 2 ⇤ k h (T h ), a h (! h , ⌘ h ) = (↵! h , ⌘ h ) ⌦ + X T 2T h (i d! h , ⌘ h ) T + (! h , j ⌘ h ) T + X f 2F @ \F @ Z f tr(i ! h ^?⌘ h ) X f 2F @ Z f tr(! h ^?j ⌘ h ) + X f 2F Z f tr({i ! h } f ^?[⌘ h ] f ) Z f tr([! h ] f ^?{j ⌘ h } f ). (16)
Proof. By using a Leibniz rule for the exterior derivative with respect to the wedge product

d(! ^?⌘) = d! ^?⌘ + ( 1) k+`! ^? ⌘ ! 2 ⇤ k (⌦) , ⌘ 2 ⇤ `(⌦)
and Stokes' theorem [START_REF] Schwarz | Hodge Decomposition-a Method for Solving Boundary Value Problems[END_REF]Theorem 1.2.7], it easily follows that

Z @⌦ tr(! ^?µ) = (d!, µ) ⌦ (!, µ) ⌦ 8! 2 ⇤ k (⌦) , µ 2 ⇤ k+1 (⌦) . (17)
Hence, using [START_REF] Lellis | Notes on hyperbolic systems of conservation laws and transport equations[END_REF] together with Cartan's homotopy formula for the adjoint of the Lie derivative L results in

X T 2T h (! h , j ⌘ h ) T = X T 2T h (i d! h , ⌘ h ) T X f 2F @ Z f tr(! h ^?j ⌘ h ) X f 2F ✓ Z f tr({! h } f ^?[j ⌘ h ] f ) + Z f tr([! h ] f ^?{j ⌘ h } f ) ◆ . ( 18
)
where the outflow boundary terms can be recast as

Z f tr(i ! h ^?⌘ h ) = Z f tr i (! h ^?⌘ h ) Z f tr(! h ^?j ⌘ h ) 8 f 2 F @ \ F @ . (19)
Therefore, substituting [START_REF] Diperna | Ordinary di↵erential equations, transport theory and Sobolev spaces[END_REF] and ( 19) into the bilinear form (13) yields the conclusion.

Note that if ⇤

k h (T h ) is a space of H⇤ k (⌦)-conforming discrete di↵erential forms, the terms tr([! h ] f ^?{j ⌘ h } f ) in (16) and cf tr([! h ] f ^?[j ⌘ h ] f ) in (14) vanish for all f 2 F and every ! h , ⌘ h 2 ⇤ k h (T h ). Remark 1. (Lipschitz continuous velocity fields 2 W 1,1 (⌦))
Let us consider the particular case of velocity fields that feature Lipschitz continuity in space, that is 2 W 1,1 (⌦). An easy computation allows to write, for all

! h , ⌘ h 2 ⇤ k h (T h ) X f 2F Z f tr({i ! h } f ^?[⌘ h ] f ) = X f 2F Z f tr(i { } f {! h } f ^?[⌘ h ] f ) + X f 2F 1 4 Z f tr(i [ ] f [! h ] f ^?[⌘ h ] f ) ; X f 2F Z f tr({! h } f ^?[j ⌘ h ] f ) = X f 2F Z f tr({! h } f ^?j { } f [⌘ h ] f ) + X f 2F Z f tr({! h } f ^?j [ ] f {⌘ h } f ) ,
and similarly for the stabilization terms in [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. Since trivially [ ] f ⌘ 0 for all f 2 F , all the terms involving the jump of the velocity can be dropped and the variational problem reduces to: find

! h 2 ⇤ k h (T h ) such that a h (! h , ⌘ h ) + s h (! h , ⌘ h ) = l(⌘ h ) for all ⌘ h 2 ⇤ k h (T h )
, where l(⌘ h ) is as in [START_REF] Brezzi | Discontinuous Galerkin methods for first-order hyperbolic problems[END_REF] while the stabilized bilinear form reads

a h (! h , ⌘ h ) + s h (! h , ⌘ h ) = (↵! h , ⌘ h ) ⌦ + X T 2T h (! h , L ⌘ h ) T 2T h + X f 2F @ \F @ Z f tr i (! h ^?⌘ h ) + X f 2F Z f tr(i {! h } f ^?[⌘ h ] f ) + Z f tr({! h } f ^?j [⌘ h ] f ) + X f 2F Z f c f tr(i [! h ] f ^?[⌘ h ] f ) + Z f cf tr([! h ] f ^?j [⌘ h ] f ).
If c f = cf , using ( 11) the bilinear form above can be recast as

a h (! h , ⌘ h ) + s h (! h , ⌘ h ) = (↵! h , ⌘ h ) ⌦ + X T 2T h (! h , L ⌘ h ) T + X f 2F @ \F @ Z f tr i (! h ^?⌘ h ) + X f 2F Z f tr i ({! h } f ^?[⌘ h ] f ) + X f 2F Z f c f tr i ([! h ] f ^?[⌘ h ] f ) (20)
and the formulation in [23, Equation 4.8, p. 61] is recovered. Note that, owing to the fact that { } f = | f , the choice of stabilization given in (15) yields a scheme with genuine upwind fluxes. Moreover, since the stabilization terms vanish for ! 2 W solution of [START_REF] Bossavit | On the geometry of electromagnetism. (2): 'Geometrical objects[END_REF], the variational formulation with stabilized bilinear form given by ( 13) and ( 14) is consistent with [START_REF] Bossavit | On the geometry of electromagnetism. (2): 'Geometrical objects[END_REF], namely Galerkin orthogonality

a h (! ! h , ⌘ h ) + s h (! ! h , ⌘ h ) = 0 8⌘ h 2 ⇤ k h (T h ) (21)
holds for ! h 2 ⇤ k h (T h ) numerical solution of the discretized problem. Observe that the stabilized Galerkin formulation [START_REF] Brezzi | Discontinuous Galerkin methods for first-order hyperbolic problems[END_REF], [START_REF] Friedrichs | Symmetric positive linear di↵erential equations[END_REF] for Lipschitz continuous velocities 2 W 1,1 (⌦) can be equivalently derived by imposing the transmission conditions (8) on the mesh facets.

Analogously, the bilinear form corresponding to the reformulated variational problem [START_REF] Davidson | An Introduction to Magnetohydrodynamics[END_REF]

for 2 W 1,1 (⌦) and ! h , ⌘ h 2 ⇤ k h (T h ) reads: a h (! h , ⌘ h ) = (↵! h , ⌘ h ) ⌦ + X T 2T h (i d! h , ⌘ h ) T + (! h , j ⌘ h ) T + X f 2F @ \F @ Z f tr(i ! h ^?⌘ h ) X f 2F @ Z f tr(! h ^?j ⌘ h ) + X f 2F Z f tr(i {! h } f ^?[⌘ h ] f ) Z f tr([! h ] f ^?j {⌘ h } f ).
2.4. Stabilized Galerkin formulation in terms of vector proxies. For the sake of completeness, we present the vector proxy representation of the stabilized reformulated bilinear form ( 16), ( 14) corresponding to the variational formulation associated with the transport problem of the corresponding k-form. Table 2 and Table 3 are used to establish the correspondences. Let V h be finite element spaces of vector proxies associated to the spaces

⇤ k h (T h ) of polynomial di↵erential k-forms on the mesh T h . Let u, v 2 V h or u, v 2 V h be the vector proxy representations of the k-forms ! h , ⌘ h 2 ⇤ k h (T h ): k = 0 : a h (u, v) + s h (u, v) = Z ⌦ ↵uv dx + X T 2T h Z T • graduv dx X f 2F @ Z f • n f uv dS + X f 2F Z f [u] f { v} f • n f dS + Z f cf [u] f [ v] f • n f dS. k = 1 : a h (u, v) + s h (u, v) = Z ⌦ ↵u • v dx + X T 2T h Z T (curl u ⇥ ) • v dx Z T u • divv dx + X f 2F @ \F @ Z f (u • )(v • n f ) dS + X f 2F @ Z f (u ⇥ n f ) • ( ⇥ v) dS + X f 2F Z f {u • } f [v] f • n f dS + Z f ([u] f ⇥ n f ) • { ⇥ v} f dS + X f 2F Z f c f [u • ] f [v] f • n f dS Z f cf ([u] f ⇥ n f ) • [ ⇥ v] f dS. k = 2 : a h (u, v) + s h (u, v) = Z ⌦ ↵u • v dx + X T 2T h Z T divu • v dx + Z T u • ( ⇥ curl v) dx + X f 2F @ \F @ Z f (u ⇥ ) • (v ⇥ n f ) dS X f 2F @ Z f (u • n f )(v • ) dS + X f 2F Z f {u ⇥ } f • ([v] f ⇥ n f ) dS Z f [u] f • n f { • v} f dS + X f 2F Z f c f [u ⇥ ] f • ([v] f ⇥ n f ) dS + Z f cf [u] f • n f [ • v] f dS. k = 3 : a h (u, v) + s h (u, v) = Z ⌦ ↵uv dx X T 2T h Z T u • gradv dx X f 2F @ \F @ Z f • n f uv dS + X f 2F Z f { u} f • n f [v] f dS + Z f c f [ u] f • n f [v] f dS.
2.5. Trial and test spaces of discrete di↵erential forms. From now, we restrict ourselves to special types of meshes:

• T h is either a simplicial decomposition of ⌦ ⇢ R n as defined in [5, Section 4.1 and Section 5.3],

• or a tensor product mesh, namely a compatible, locally quasi-uniform, a ne mesh partition of ⌦ into non-degenerate axiparallel parallelotopes.

On such meshes various families of piecewise polynomial discrete di↵erential forms of any degree have been constructed, see [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF], [START_REF] Arnold | Spaces of finite element di↵erential forms[END_REF], [27, Section 3] and [START_REF] Arnold | Finite element di↵erential forms on curvilinear cubic meshes and their approximation properties[END_REF] for a detailed overview.

For a simplicial decomposition T h , the spaces of polynomial totally discontinuous discrete di↵erential k-forms on T h are defined as

P d r ⇤ k (T h ) := {! 2 L 2 ⇤ k (⌦) , ! | T 2 P r ⇤ k (T ), T 2 T h }
where P r ⇤ k (T ) is the space of di↵erential k-forms with polynomial coe cients of degree at most r on the n-cell T 2 T h , obtained as the restriction of

P r ⇤ k (R n ) to T . The corresponding space of H⇤ k (⌦)-conforming discrete di↵erential forms P r ⇤ k (T h ) := {! 2 H⇤ k (⌦) , ! | T 2 P r ⇤ k (T ), T 2 T h } allows to introduce another family of polynomial di↵erential k-forms on T h , namely P r ⇤ k (T h ) = P r 1 ⇤ k (T h ) H r 1 ⇤ k+1 (T h ) = P r 1 ⇤ k (T h ) + P r 1 ⇤ k+1 (T h ) or equivalently P r ⇤ k (T h ) := {! 2 P r ⇤ k (T h ) : ! 2 P r ⇤ k 1 (T h )}, where H r ⇤ k (T h )
is the space of homogeneous polynomial di↵erential k-forms of degree r and  :

H r ⇤ k (T h ) ! H r+1 ⇤ k 1 (T h ) denotes the Koszul di↵erential [5, Section 3.2].
The so-called "first family" of finite element di↵erential k-forms is hence defined as

P r ⇤ k (T h ) = {! 2 H⇤ k (⌦) : ! | T 2 P r ⇤ k (T )
, T 2 T h } whose functions satisfy the continuity requirement that the trace tr ! is single-valued on all n 1-cells which in turn ensures inclusion in H⇤ k (⌦). The family Q r ⇤ k (T h ) of finite element di↵erential forms on a tensor product mesh can be constructed by iteratively applying a tensor product strategy from the 1-dimensional interval. This tensor product construction allows to build a subcomplex of the de Rham complex. We refer the interested reader to [4, Section 5] for details on such constructions.

The H⇤ k (⌦)-conforming finite element spaces presented above over a cell complex T h form a discrete de Rham sequence as a cochain projection from the de Rham complex through projection operators ⇧ k h , namely the following diagram commutes

H⇤ 0 (⌦) d ! H⇤ 1 (⌦) d ! . . . d ! H⇤ n (⌦) ? ? y⇧ 0 h ? ? y⇧ 1 h ? ? y ⇧ n h ⇤ 0 h (T h ) d ! ⇤ 1 h (T h ) d ! . . . d ! ⇤ n h (T h ) (22)
where each [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF]Section 5.5] and [START_REF] Arnold | Spaces of finite element di↵erential forms[END_REF]).

⇤ k h (T h ) ! ⇤ k+1 h (T h ) can be substituted by P r ⇤ k (T h ) ! P r 1 ⇤ k+1 (T h ), or P r ⇤ k (T h ) ! P r ⇤ k+1 (T h ), or Q r ⇤ k (T h ) ! Q r ⇤ k+1 (T h ) for every r 1 (see
In the subsequent analysis we will make use of pairs of H⇤ k (⌦)-conforming spaces ⇤ k h (T h ) and non-conforming spaces ⇤ d,k h (T h ) as in the following:

(I) ⇤ d,k h (T h ) = P d r ⇤ k (T h ) and ⇤ k h (T h ) = P r ⇤ k (T h ) with T h simplicial mesh; (II) ⇤ d,k h (T h ) = {! 2 L 2 ⇤ k (⌦) , ! | T 2 P r+1 ⇤ k (T ), T 2 T h } and ⇤ k h (T h ) = P r+1 ⇤ k (T h ) with T h simplicial mesh; (III) ⇤ d,k h (T h ) = {! 2 L 2 ⇤ k (⌦) , ! | T 2 Q r+1 ⇤ k (T ), T 2 T h } and ⇤ k h (T h ) = Q r+1 ⇤ k (T h ) with T h tensor product mesh.
Remark 2. We pay particular attention to H⇤ k (⌦)-conforming trial/test spaces because they allow the straightforward Galerkin discretization of the di↵usion form d ? d present in (3).

3. Stationary transport: Estimates for continuous velocity fields. As explained in Section 1.2, a rigorous convergence analysis is only possible in the case 2 W 1,1 (⌦), for want of a well-posedness result for [START_REF] Bossavit | On the geometry of electromagnetism. (2): 'Geometrical objects[END_REF] with discontinuous velocity fields. Hence, all theoretical results will rely on the assumption 2 W 1,1 (⌦). Moreover, without loss of generality we can assume the scaling

k k L 1 (⌦) = 1.
Also, in this section, we admit only the same special classes of meshes as in Section 2.5. We recall the notion of shape regularity of a mesh as presented, e.g., in [START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF]Section 3.6] following [14, Section 3.1]. Throughout, let ⇤ k h (T h ) denote some piecewise polynomial approximation space for di↵erential k-forms. If not otherwise specified, ⇤ k h (T h ) could be either a H⇤ k (⌦)-conforming approximation space P r ⇤ k (T h ) or P r+1 ⇤ k (T h ), but also the totally discontinuous space P d r ⇤ k (T h ) on a simplicial mesh and Q r+1 ⇤ k (T h ) on a tensor product mesh. Note that we will always assume that [START_REF] Friedrichs | Symmetric positive linear di↵erential equations[END_REF]. Note that the bilinear form s h (•, •) associated to the stabilization operator is symmetric and non-negative on V (h)⇥V (h). Moreover, for all ⌘ h 2 ⇤ k h (T h ), applying (9) to a h (⌘ h , ⌘ h ), results in

⇤ n h (T h ) = P d r ⇤ n (T h ). Let V (h) := ⇤ k h (T h ) + V , we introduce the discrete operators A h , S h : V (h) ! ⇤ k h (T h ) such that for all ! 2 V (h), ⌘ h 2 ⇤ k h (T h ), (A h !, ⌘ h ) ⌦ := a h (!, ⌘ h ) and (S h !, ⌘ h ) ⌦ := s h (!, ⌘ h ) , with a h (•, •) and s h (•, •) as in
a h (⌘ h , ⌘ h ) = (↵⌘ h , ⌘ h ) ⌦ + X T 2T h (⌘ h , L ⌘ h ) T + X f 2F @ \F @ Z f tr i (⌘ h ^?⌘ h ) + X f 2F Z f tr i ({⌘ h } f ^?[⌘ h ] f ) = X T 2T h ✓ ⌘ h , ↵ + 1 2 (L + L )⌘ h ◆ T + 1 2 X f 2F \F @ Z f tr i (⌘ h ^?⌘ h ) 1 2 X f 2F @ Z f tr i (⌘ h ^?⌘ h ) = 1 2 Z @⌦ | • n @⌦ | tr i n @⌦ (⌘ h ^?⌘ h ) 1 2 (⇤⌘ h , ⌘ h ) ⌦ (23)
where ⇤ := (2↵ id +L + L ). Let us introduce the following norms on V (h),

k!k 2 h := k!k 2 L 2 ⇤ k (⌦) + |!| 2 h ; with |!| 2 h := X f 2F @ \F @ k!k 2 f, + X f 2F @ k!k 2 f, + X f 2F k[!] f k 2 f,c f , (24)
k!k 2 f, := Z f tr i (! ^?!) , and k!k 2 f,c f := Z f c f tr i (! ^?!) .
Note that the above norms are well-defined in view of the definition of inflow and outflow boundary and the fact that tr i (!

^?!) = ( • n f ) tr i n f (! ^?!) together
with Assumption 2 on the stabilization coe cients c f , f 2 F . The above norms are combined into

k!k ⇤ := h 1/2 k!k L 2 ⇤ k (⌦) + h 1/2 |!| H 1 ⇤ k (T h ) + ✓ X f 2F k!k 2 L 2 ⇤ k (f ) ◆ 1/2 + |!| h ( 25
)
where

|•| H 1 ⇤ k (T h
) stands for the broken H⇤ k -seminorm on T h . Convergence estimates for the spatial discretization of the stationary boundary value problem are key to analyzing the convergence of the fully discrete scheme. They hinge on stability results for the di↵erential operator L h := A h + S h . In order for these results to hold for both non-conforming and H⇤ k (⌦)-conforming space discretization, we approximate discontinuous di↵erential forms by di↵erential forms in H⇤ k (T h ) as in the following:

Proposition 2. Let ⇤ d,k h (T h
) and ⇤ k h (T h ) be defined as in either (I), (II) or (III).

Then for every

! 2 ⇤ d,k h (T h ), there exists ! c 2 ⇤ k h (T h ) such that k! ! c k 2 L 2 ⇤ k (⌦)  Ch X f 2F ktr[!]k 2 L 2 ⇤ k (f )
with constant C > 0 depending only on the polynomial degree and the shape regularity of the mesh.

The proof of this result is constructive and involves some technicalities. We omit it here and refer the interested reader to [26, Appendix A]. Note that the construction of the H⇤ k (⌦)-conforming approximation is based on an averaging interpolation which is the extension to discrete di↵erential k-forms of the operator introduced and studied for scalar functions in R 3 in [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems[END_REF] and [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF]Appendix]. , the polynomial degree r and the shape regularity of the mesh, such that

|! h | h  C S h 1/2 k! h k L 2 ⇤ k (⌦) 8 ! h 2 ⇤ k h (T h ) . (26)
Moreover, for every

! 2 V (h) kL h !k L 2 ⇤ k (⌦)  C L k!k L 2 ⇤ k (⌦) + |!| H 1 ⇤ k (T h ) + C 0 L h 1/2 |!| h ; ( 27
)
where the constant C L depends on ↵ and | | W 1,1 (⌦) , and the constant C 0 L depends on the stabilization coe cients

|c f | 1/2 , |c f | 1/2
, the polynomial degree r and the shape regularity of the mesh.

Proof. The first inequality [START_REF] Heumann | Stabilized Galerkin for Transient Advection of Di↵erential Forms[END_REF] immediately follows by the definition of h-seminorm in [START_REF] Heumann | Stabilized Galerkin methods for magnetic advection[END_REF] and inverse trace inequalities [14, p. 146]. The proof of ( 27) is based on standard norm inequalities and can be found in [START_REF] Heumann | Stabilized Galerkin for Transient Advection of Di↵erential Forms[END_REF]Lemma 3.2].

For the sake of conciseness, the following bounds on orthogonal subscales are presented in the case of full polynomial spaces of discrete di↵erential k-forms on simplices, case (I), namely for the spaces

⇤ d,k h (T h ) = P d r ⇤ k (T h ) and ⇤ k h (T h ) = P r ⇤ k (T h ).
The proof for the cases (II) and on tensor product meshes (III), follows mutatis mutandis.

Lemma 3.2. The following statements hold true:

(i) If ⇧ h denotes the L 2 -orthogonal projection onto ⇤ d,k h (T h ), then there exists a constant C ⇡ such that for all ! 2 ⇤ d,k h (T h ) + V , ⌘ h 2 ⇤ d,k h (T h ) with r 1, it holds |(L h (! ⇧ h !), ⌘ h ) ⌦ |  C ⇡ k! ⇧ h !k ⇤ k⌘ h k h ; (ii) If ⇧ h denotes the global L 2 -orthogonal projection onto ⇤ k h (T h ), then there exists a constant C 0 ⇡ such that for all ! 2 V (h), ⌘ h 2 ⇤ k h (T h ) |(L h (! ⇧ h !), ⌘ h ) ⌦ |  C 0 ⇡ k! ⇧ h !k ⇤ k⌘ h k h ;
where the constants C ⇡ and C 0 ⇡ depend on ↵, | | W 1,1 (⌦) , the stabilization coef-

ficients |c f | 1/2 , |c f | 1/2
, the polynomial degree r and the shape regularity of the mesh.

Proof. In order to show (i), we proceed as in [24, Theorem 3.1]. Let h be the L 2 -projection of 2 W 1,1 (⌦) onto piecewise constant vector fields. One can add the zero term [START_REF] Friedrichs | Symmetric positive linear di↵erential equations[END_REF]. Hence, using estimates on the projection error for in the L 1 -norm, Cauchy-Schwarz inequality and inverse inequalities, results in

P T 2T h ! ⇧ h !, L h ⌘ h T to the bilinear form (L h (! ⇧ h !), ⌘ h ) ⌦ given in
|(L h (! ⇧ h !), ⌘ h ) ⌦ |  k↵k L 1 (⌦) k! ⇧ h !k L 2 ⇤ k (⌦) k⌘ h k L 2 ⇤ k (⌦) + | h | W 1,1 (⌦) k! ⇧ h !k L 2 ⇤ k (⌦) k⌘ h k L 2 ⇤ k (⌦) + C| | W 1,1 (⌦) k! ⇧ h !k L 2 ⇤ k (⌦) k⌘ h k L 2 ⇤ k (⌦) + X f 2F k[⌘ h ]k f,c f k[! ⇧ h !] c 1 f {! ⇧ h !}k f,c f + |⌘ h | h X f 2F @ \F @ k! ⇧ h !k f,  C ⇡ k! ⇧ h !k ⇤ k⌘ h k h 8 ⌘ h 2 ⇤ d,k h (T h ) , (28)
where the interior facet terms above have been bounded as follows. Let µ := ! ⇧ h !. Using inverse trace inequalities results in

X f 2F k[µ] c 1 f {µ}k 2 f,c f  X f 2F f =@T + \@T |c f |k k L 1 (⌦) ✓ kµk 2 L 2 ⇤ k (@T + ) + kµk 2 L 2 ⇤ k (@T ) + kc 1 f µk 2 L 2 ⇤ k (@T + ) + kc 1 f µk 2 L 2 ⇤ k (@T ) ◆  X f 2F f =@T + \@T |c f |k k L 1 (⌦) max{1, |c f | 2 }kµk 2 L 2 ⇤ k (@T + [@T )  k k L 1 (⌦) max f 2F max{|c f |, |c f | 1 } X T 2T h kµk 2 L 2 ⇤ k (@T ) .
In the case (ii) of H⇤ k (⌦)-conforming discretization, we can proceed as above, use estimate [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF], but the non-zero term

P T 2T h ! ⇧ h !, L h ⌘ h T has to be bounded. We show that for all ! 2 V (h), ⌘ h 2 P r ⇤ k (T h ), X T 2T h ! ⇧ h !, L h ⌘ h T  Ch 1/2 k! ⇧ h !k L 2 ⇤ k (⌦) k⌘ h k h , (29)
with the constant C > 0 depending only on the polynomial degree and the shape regularity of the mesh. In order to do that, we exploit the fact that since h is piecewise constant for all T 2 T h , L h = L h and we build H⇤ k (⌦)-conforming approximations for each of the two terms appearing in Cartan's formula L h = i h d + di h . In particular, in view of Proposition 2, let c,k 2 P r ⇤ k (T h ) be the H⇤ k (⌦)-conforming approximation of i h d⌘ h 2 P d r ⇤ k (T h ) and let c,k 1 2

P r+1 ⇤ k 1 (T h ) be the H⇤ k (⌦)-conforming approximation of i h ⌘ h 2 P d r+1 ⇤ k 1 (T h ). Since c,k , d c,k 1 2 P r ⇤ k (T h ), ! ⇧ h !, L h ⌘ h T = ! ⇧ h !, i h d⌘ h c,k T + ! ⇧ h !, d(i h ⌘ h c,k 1 ) T
and by Cauchy-Schwarz inequality and inverse inequality, one has X

T 2T h ! ⇧ h !, L h ⌘ h T  k! ⇧ h !k L 2 ⇤ k (⌦) i h d⌘ h c,k L 2 ⇤ k (⌦) + Ch 1 k! ⇧ h !k L 2 ⇤ k (⌦) i h ⌘ h c,k 1 L 2 ⇤ k 1 (⌦)
. By the approximation results in Proposition 2, the projection errors can be bounded as X

T 2T h ! ⇧ h !, L h ⌘ h T  Ch k! ⇧ h !k L 2 ⇤ k (⌦) ✓ X f 2F ktr ⇥ i h d⌘ h ⇤ f k 2 L 2 ⇤ k (f ) ◆ 1/2 + k! ⇧ h !k L 2 ⇤ k (⌦) ✓ X f 2F ktr ⇥ i h ⌘ h ⇤ f k 2 L 2 ⇤ k 1 (f ) ◆ 1/2 .
Upper bounds for the facet terms can be derived as follows: Let us decompose the velocity field into its normal component n := ( • n)n and its tangential component t := (n ⇥ ) ⇥ n. Then we can write

tr(i ⌘ h ) = tr(i n ⌘ h + i t ⌘ h ) = ( • n) tr(i n ⌘ h ) + i t tr ⌘ h 8⌘ h 2 ⇤ k h (T h ) , 8k. (30
) If f = @T + \ @T , using estimates on the projection error for , trace and inverse inequalities together with [START_REF] Kato | Linear and Quasi-Linear Equations of Evolution of Hyperbolic Type, Hyperbolicity[END_REF] and the fact that d⌘ h 2 H⇤ k+1 (⌦) due to [START_REF] Goedbloed | Principles of Magnetohydrodynamics[END_REF], results in

ktr ⇥ i h d⌘ h ⇤ f k L 2 ⇤ k (f )  ktr ⇥ i ( h ) d⌘ h ⇤ f k L 2 ⇤ k (f ) + ktr [i d⌘ h ] f k L 2 ⇤ k (f )  k h k L 1 (⌦) k[d⌘ h ] f k L 2 ⇤ k+1 (f ) + k( • n f ) tr i n f [d⌘ h ] f k L 2 ⇤ k (f )  Ch| | W 1,1 (⌦) h 1/2 h 1 k⌘ h k L 2 ⇤ k (T + [T ) + Ch 1 k[⌘ h ] f k f, . Similarly, ⌘ h 2 P r ⇤ k (T h ) ⇢ H⇤ k (⌦) implies ⇥ i t tr ⌘ h ⇤ f = 0 for all f 2 F, hence ktr ⇥ i h ⌘ h ⇤ f k L 2 ⇤ k 1 (f )  ktr ⇥ i ( h ) ⌘ h ⇤ f k L 2 ⇤ k 1 (f ) + ktr [i ⌘ h ] f k L 2 ⇤ k 1 (f )  Ch| | W 1,1 (⌦) h 1/2 k⌘ h k L 2 ⇤ k (T + [T ) + k[⌘ h ] f k f,
which leads to the desired estimate [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems[END_REF]. Finally, combining the estimates ( 28) and (29) yields

| (L h (! ⇧ h !), ⌘ h ) ⌦ |  C ⇡ k! ⇧ h !k ⇤ k⌘ h k h + Ch 1/2 k! ⇧ h !k L 2 ⇤ k (⌦) k⌘ h k h  C 0 ⇡ k! ⇧ h !k ⇤ k⌘ h k h .
Note that [START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF] together with the definition of k•k ⇤ in [START_REF] Heumann | Convergence of Lowest Order Semi-Lagrangian Schemes[END_REF], inverse and trace inequalities gives

kL h !k L 2 ⇤ k (⌦)  Ch 1/2 k!k ⇤  Ch 1 k!k L 2 ⇤ k (⌦) 8! 2 V (h), (31)
and C depends only on the polynomial degree and the shape regularity of the mesh. As it will be shown in the following, stability of second order Runge-Kutta schemes can be achieved with the standard CFL condition if the space discretization is performed with piecewise linear finite elements. Therefore, this case is tackled separately. In particular, we can establish the following estimate.

Lemma 3.3. Let ⇧ 0 h denote the L 2 -orthogonal projection onto P d 0 ⇤ k (T h ).
In the case of space discretization with piecewise a ne elements

⇤ k h (T h ) = P d 1 ⇤ k (T h ), or ⇤ k h (T h ) = P 1 ⇤ k (T h ) or ⇤ k h (T h ) = P 1 ⇤ k (T h ) or ⇤ k h (T h ) = Q 1 ⇤ k (T h ), there exists a constant C ⇡ which depends on ↵, | | W 1,1 (⌦) , the stabilization coe cients |c f | 1/2 , |c f | 1/2
and the shape regularity of the mesh, such that for all

! h , ⌘ h 2 ⇤ k h (T h ) |(L h ! h , ⌘ h ⇧ 0 h ⌘ h ) ⌦ |  C ⇡ h 1/2 k! h k h ⌘ h ⇧ 0 h ⌘ h L 2 ⇤ k (⌦) .
Proof. The idea of the proof is very similar to the one proposed in [13, Lemma 2.1] and follows a reasoning analogous to the one in Theorem 3.

a h (!, !) + s h (!, !) min ⇢ 1 2 ↵ 0 , 1 k!k 2 h 8! 2 ⇤ k h (T h ) . ( 32
)
Moreover, if ! 2 H r+1 ⇤ k (⌦) is solution of the advection problem (5) and ! h 2 ⇤ k h (T h ) is solution of the discrete variational formulation with bilinear form given in (20), then k! ! h k h  Ch r+1/2 k!k H r+1 ⇤ k (⌦)
with the constant C > 0 depending on |c f |, |c f | 1 , ↵, , the polynomial degree r and the shape regularity of the mesh.

Proof. The proof of stability [START_REF] Levy | From semidiscrete to fully discrete: stability of Runge-Kutta schemes by the energy method[END_REF] immediately follows by [START_REF] Heumann | Eulerian and Semi-Lagrangian Methods for Advection-Di↵usion of Di↵erential Forms[END_REF], the positivity condition [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF] and the definition of the h-norm.

Let ⇧ h denote the L 2 -projection into ⇤ k h (T h ). By stability and consistency ( 21), one has min

⇢ 1 2 ↵ 0 , 1 k! ⇧ h !k 2 h  |(L h (! ⇧ h !), ⌘ h ) ⌦ |
where ⌘ h := ! h ⇧ h !. We proceed as in the proof of Theorem 3.2 (i) to get [START_REF] Houston | Interior penalty method for the indefinite time-harmonic Maxwell equations[END_REF] and use a multiplicative trace inequality (see [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]Theorem 1.6.6]) for the interior facets terms, i.e.

k! ⇧ h !k 2 L 2 ⇤ k (@T )  C ⇣ h 1 T k! ⇧ h !k 2 L 2 ⇤ k (T ) + h T |! ⇧ h !| 2 H 1 ⇤ k (T )
⌘ with C depending only on the shape of T . Moreover, in the case of H⇤ k (⌦)conforming discrete di↵erential forms the extra non-zero terms are bounded as in [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems[END_REF]. The approximation estimates [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF]Theorem 5.3] inf

µ h 2Pr⇤ k (T ) k! µ h k L 2 ⇤ k (T )  Ch r+1 k!k H r+1 ⇤ k (T ) ; inf µ h 2Pr⇤ k (T ) k! µ h k H 1 ⇤ k (T )  Ch r k!k H r+1 ⇤ k (T ) ;
for C > 0 independent of h, yield the conclusion.

Note that in the non-stabilized case (c f = 0), the bilinear form in ( 20) is coercive in the L 2 ⇤ k -norm but only sub-optimal convergence is attained, namely

k! ! h k L 2 ⇤ k (⌦)  Ch r k!k H r+1 ⇤ k (⌦)
holds with C > 0 independent of the mesh width h.

4.

Fully discrete problem. In the present section, we formulate the fully discrete advection problem for a di↵erential k-form by coupling the stabilized Galerkin spatial discretization introduced in Section 2 with explicit timestepping schemes. In particular, the forward Euler method and explicit second-order and third-order Runge-Kutta (RK) schemes are investigated.

On the time interval I = [0, T ], we consider a uniform partition S N 1 n=0 [t n , t n+1 ] for a given positive integer N and t n = n⌧ with uniform time step ⌧ such that T = N⌧ . The semi-discrete problem reads: find

! h (t) 2 ⇤ k h (T h ) such that (@ t ! h (t), ⌘ h ) ⌦ + (L h (t)! h (t), ⌘ h ) ⌦ = l(t)(⌘ h ) 8⌘ h 2 ⇤ k h (T h ) (! h (0), ⌘ h ) ⌦ = ! 0 , ⌘ h ⌦ 8⌘ h 2 ⇤ k h (T h ) (33)
where the bilinear forms a h (•, •), s h (•, •) and the linear functional l(•) are obtained at each time step through spatial discretization as in ( 13), ( 14) and ( 12) with forcing term f (x, t) and velocity field (x, t). The semi-discrete problem [START_REF] Mishra | Multi-level Monte Carlo finite volume methods for nonlinear systems of conservation laws in multi-dimensions[END_REF] can equivalently be recast as the finite dimensional operator evolution equation

@ t ! h (t) + L h (t)! h (t) = F h (t) 8t 2 [0, T ] ( 34
)
where

F h (t) 2 ⇤ k h (T h ) is such that (F h (t), ⌘ h ) ⌦ = l(t)(⌘ h ) for all ⌘ h 2 ⇤ k h (T h ).
In light of the results established in Theorem 3.1, Theorem 3.2 and Theorem 3.3, quasi-optimal convergence rates for the L 2 ⇤ k -error in space L 1 -error in time can be proven for smooth solutions of the problem [START_REF] Mishra | Multi-level Monte Carlo finite volume methods for nonlinear systems of conservation laws in multi-dimensions[END_REF], along the lines of the analysis proposed by Burman et al. in [START_REF] Burman | Explicit Runge-Kutta schemes and finite elements with symmetric stabilization for first-order linear PDE systems[END_REF]. In particular, under CFL-type conditions, the e cacy of the proposed space-time discretization lies in the fact that the anti-di↵usive nature of explicit RK schemes is compensated by the artificial dissipation introduced through the stabilized spatial discretization.

In the following paragraphs, we introduce the fully discrete problem for (34) by explicitly stating the stages corresponding to the time-stepping. Moreover, we present the convergence results corresponding to each fully discrete scheme. 4.1. Explicit Euler Scheme. The first order explicit Euler scheme for the problem [START_REF] Orszag | Small-scale structure of two-dimensional magnetohydrodynamic turbulence[END_REF] reads

! n+1 h = ! n h ⌧ L n h ! n h + ⌧ F n h ( 35
)
where

! n h = ! h (t n ), L n h := L h (t n ) and F n h := F h (t n ). Theorem 4.1. Let ! 2 C 0 (0, T ; H r+1 ⇤ k (⌦)) \ C 2 (0, T ; L 2 ⇤ k (⌦)) be the exact solution of (5) and let {! n h } N n=1 ⇢ ⇤ k h (T h
) be the discrete solution of problem [START_REF] Pazy | Semigroups of linear operators and applications to partial di↵erential equations[END_REF]. Let Assumption 2 and the monotonicity condition (6) for 2 C 0 (0, T ; W 1,1 (⌦)) hold true. Consider the trial spaces

⇤ k h (T h ) = P d r ⇤ k (T h ) or ⇤ k h (T h ) = P r ⇤ k (T h ) or ⇤ k h (T h ) = P r+1 ⇤ k (T h ) or ⇤ k h (T h ) = Q r+1 ⇤ k (T h ).
Then there exist constants C, C CFL > 0 depending only on the constants in Theorem 3.1, Theorem 3.2 and Theorem 3.3 and the trial/test spaces

⇤ k h (T h ) such that, if (i) ⌧  C CFL h, for ⇤ k h (T h ) = P d 0 ⇤ k (T h ); (ii) ⌧  C CFL h 2
, for all other choices of ⇤ k h (T h );

then max 0nN k!(t n ) ! n h k L 2 ⇤ k (⌦)  C(⌧ + h r+1/2 ).
Proof. For the sake of clarity we briefly sketch here the underlying idea proposed in [START_REF] Levy | From semidiscrete to fully discrete: stability of Runge-Kutta schemes by the energy method[END_REF] and [START_REF] Burman | Explicit Runge-Kutta schemes and finite elements with symmetric stabilization for first-order linear PDE systems[END_REF], valid also for the higher order timestepping schemes introduced below. For ! n = !(t n ) exact solution at the n-th time step, the proof starts with writing the error generated at each stage (here one stage) of the scheme as

! n ! n h = (! n ⇧ h ! n ) (! n h ⇧ h ! n ) =: e n ⇡ e n h
and bounding the error e n h by the approximation error e n ⇡ . In order to do that, by deriving the equation governing the time evolution of the error e n h , an energy identity associated with the timestepping scheme can be identified. Starting from this identity, the desired estimate is obtained by deriving upper bounds on each terms via the estimates established in Theorem 3.1, Theorem 3.2 and Theorem 3.3. Note that [START_REF] Lee | Introduction to smooth manifolds[END_REF] is crucial in order to achieve stability under the CFL-condition in (ii), while the bound on the orthogonal subscales inferred in Theorem 3.2 and Theorem 3.3 is instrumental in the derivation of quasi-optimal error estimates. The conclusion follows by a Gronwall type argument and standard estimates on the projection error e n ⇡ .

Remark 3. Note that the mild time step constraint in Theorem 4.1 (i) valid for spatial approximations based on piecewise constants discontinuous elements (finite volume)

⇤ k h (T h ) = P d 0 ⇤ k (T h ) hinges on the trivial observation that the bound (27) in Theorem 3.1 reduces to kL h !k L 2 ⇤ k (⌦)  C L k!k L 2 ⇤ k (⌦) + C 0 L h 1/2 |!| h for all ! 2 V (h)
. This is the standard CFL-condition for upwind finite volume or finite di↵erence schemes for scalar advection problems. 4.2. Explicit RK2 Schemes. We consider, as in [13, Section 3.1], explicit Runge-Kutta scheme of order two (RK2) for the problem (34) of the form

µ n h = ! n h ⌧ L n h ! n h + ⌧ F n h (36)
! n+1 h = 1 2 (µ n h + ! n h ) 1 2 ⌧ L n h µ n h + 1 2 ⌧ n h (37)
where

n h := F n h + ⌧ (@ t F h )(t n ) + n h , for f in (5) su ciently smooth in time and n h such that k n h k L 2 ⇤ k (⌦)  C⌧ 2 .
Similarly to the explicit Euler scheme, convergence of the fully discrete problem with second order two-stage Runge-Kutta schemes of the form ( 36), ( 37) can be established. The proof of the following theorem is similar to that of [13, Theorem 3.1 and Theorem 3.2].

Theorem 4.2. Let ! 2 C 0 (0, T ; H r+1 ⇤ k (⌦)) \ C 3 (0, T ; L 2 ⇤ k (⌦)) be the ex- act solution of (5) and let {! n h } N n=1 ⇢ ⇤ k h (T h
) be the discrete solution of problem (36)- [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF]. Let Assumption 2 and the monotonicity condition (6) for 2 C 0 (0, T ; W 1,1 (⌦)) hold true. Consider the trial spaces

⇤ k h (T h ) = P d r ⇤ k (T h ) or ⇤ k h (T h ) = P r ⇤ k (T h ) or ⇤ k h (T h ) = P r+1 ⇤ k (T h ) or ⇤ k h (T h ) = Q r+1 ⇤ k (T h ).
Then there exist constants C, C CFL > 0 depending only on the constants in Theorem 3.1, Theorem 3.2 and Theorem 3.3 and the trial/test spaces

⇤ k h (T h ) such that, if (i) ⌧  C CFL h, for a non-conforming discretization with ⇤ k h (T h ) = P d 0 ⇤ k (T h ) or ⇤ k h (T h ) = P d 1 ⇤ k (T h ) or for a H⇤ k (⌦)-conforming approximation with spaces ⇤ k h (T h ) = P 1 ⇤ k (T h ) or ⇤ k h (T h ) = P 1 ⇤ k (T h ) or ⇤ k h (T h ) = Q 1 ⇤ k (T h ); (ii) ⌧  C CFL h 4/3
, for all other choices of ⇤ k h (T h );

then max 0nN k!(t n ) ! n h k L 2 ⇤ k (⌦)  C(⌧ 2 + h r+1/2 ).
4.3. Explicit RK3 Schemes. The explicit Runge-Kutta scheme of order three (RK3) for the problem [START_REF] Orszag | Small-scale structure of two-dimensional magnetohydrodynamic turbulence[END_REF] as in [13, Section 4.1] reads

µ n h = ! n h ⌧ L n h ! n h + ⌧ F n h (38)
n h = 1 2 (µ n h + ! n h ) 1 2 ⌧ L n h µ n h + 1 2 ⌧ (F n h + ⌧ (@ t F h )(t n )) (39)
! n+1 h = 1 3 (µ n h + n h + ! n h ) 1 3 ⌧ L n h µ n h + 1 3 ⌧ n h ( 40
)
where

n h := F n h +⌧ (@ t F h )(t n )+ 1 2 ⌧ 2 (@ tt F h )(t n )+ n h with f in (5) su ciently smooth in time and k n h k L 2 ⇤ k (⌦)  C⌧ 3 .
The proof of the following theorem follows that of [START_REF] Burman | Explicit Runge-Kutta schemes and finite elements with symmetric stabilization for first-order linear PDE systems[END_REF]Theorem 4.1].

Theorem 4.3. Let ! 2 C 0 (0, T ; H r+1 ⇤ k (⌦)) \ C 4 (0, T ; L 2 ⇤ k (⌦)) be the exact solution of (5) and let {! n h } N n=1 ⇢ ⇤ k h (T h ) be the discrete solution of problem (38)
, ( 39) and (40). Let Assumption 2 and the monotonicity condition (6) for 2 C 0 (0, T ; W 1,1 (⌦)) hold true. Consider the trial spaces

⇤ k h (T h ) = P d r ⇤ k (T h ) or ⇤ k h (T h ) = P r ⇤ k (T h ) or ⇤ k h (T h ) = P r+1 ⇤ k (T h ) or ⇤ k h (T h ) = Q r+1 ⇤ k (T h ).
Then there exist constants C, C CFL > 0 depending only on the constants in Theorem 3.1, Theorem 3.2 and Theorem 3.3 and the trial/test spaces

⇤ k h (T h ) such that, under the 1-CFL condition ⌧  C CFL h it holds max 0nN k!(t n ) ! n h k L 2 ⇤ k (⌦)  C(⌧ 3 + h r+1/2 ).

5.

Numerical experiments in 2D: continuous velocity. The two dimensional transport problem written for a 1-form amounts of solving (as in [24, Section 5])

@ t u + ↵u + grad( • u) R div(Ru) = f in ⌦ ⇥ [0, T ] u = g on in ⇥ [0, T ] u(0) = u 0 in ⌦ (41)
where R corresponds to a clockwise rotation of ⇡/2. We consider the pure transport (↵ = 0) problem in the domain ⌦ = [ 1, 1] 2 and in the time interval [0, 1]. The time-dependent Lipschitz continuous velocity field is

(x, t) := sin(2⇡t)(1 x 2 )(1 y 2 ) ✓ y x
◆ so that there is no inflow boundary. The initial condition is defined as the "bump"

u 0 := ⇢ (', ') > if x 2 + (y 0.25) 2 < 0.25 (0, 0) > otherwise (42) with '(x, y) := ✓ cos(⇡ p x 2 + (y 0.25) 2 ) 4 cos(⇡ p x 2 + (y 0.25) 2 ) 4 ◆
and the forcing term is set to zero, f = (0, 0) > . Since at time t = 1 the exact solution coincides with the initial condition u 0 in (42), we compare the numerical solution obtained at the final time T = 1 with u 0 . Figure 1 shows the L 2 -error at the final time, for a numerical spatial discretization based on non-conforming fully discontinuous piecewise polynomial vector-valued functions (a) and H(curl, ⌦)-conforming rotated Raviart-Thomas [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF] elements (b) of polynomial degree r and explicit Euler timestepping. The CFL-condition is chosen in order to fulfill the hypotheses of Theorem 4.1. Note that the error of the fully discrete scheme attains the convergence rates predicted by the theory. Analogously, the convergence rates reported for the Heun method (Figure 2) and for the RK3 scheme (Figure 3) comply with the error behavior derived in Theorem 4.2 and Theorem 4.3, respectively. Remark 4. Though we observed an increased convergence rate of the spatial discretization (also in the stationary case) compared to the predictions of the theory, our results are probably sharp. In the case of scalar transport it is well-known [START_REF] Peterson | A note on the convergence of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF] that on very special meshes, sometimes called Peterson meshes, the L 2 -norm estimates are sharp. -error -error -error

10 1 h L 2 -error r=0 ⌧ =O(h) r=1 ⌧ =O(h 2 ) O(⌧ + h r+1 ) (a) ⇤ 1 h (T h ) = P d r ⇤ 1 (T h ) 10 2 10 1 10 3 10 2 10 1 h L 2 -error r=0 ⌧ =O(h 2 ) r=1 ⌧ =O(h 2 ) O(⌧ + h r+1 ) (b) ⇤ 1 h (T h ) = P r+1 ⇤ 1 (T h )
r=0 ⌧ =O(h) r=1 ⌧ =O(h) r=2 ⌧ =O(h 4/3 ) r=3 ⌧ =O(h 4/3 ) O(⌧ 2 + h r+1 ) (a) ⇤ 1 h (T h ) = P d r ⇤ 1 (T h )
r=0 ⌧ =O(h) r=1 ⌧ =O(h 4/3 ) r=2 ⌧ =O(h 4/3 ) O(⌧ 2 + h r+1 ) (b) ⇤ 1 h (T h ) = P r+1 ⇤ 1 (T h )
r=0 ⌧ =O(h) r=1 ⌧ =O(h) r=2 ⌧ =O(h) O(⌧ 3 + h r+1 ) (a) ⇤ 1 h (T h ) = P d r ⇤ 1 (T h ) 10 2 10 1 10 6 10 4 10 2 h L 2 -error r=0 ⌧ =O(h) r=1 ⌧ =O(h) r=2 ⌧ =O(h) O(⌧ 3 + h r+1 ) (b) ⇤ 1 h (T h ) = P r+1 ⇤ 1 (T h )
Figure 3. Stabilized Galerkin discretization with upwind stabilization, polynomial di↵erential 1-forms ⇤ 1 h (T h ) and third order explicit Runge-Kutta timestepping.

6.

Numerical experiments in 2D: discontinuous velocity. Lacking a sound convergence theory for the generalized advection problem in the presence of discontinuous velocity fields, the first set of experiments aims at testing the numerical performances of the stabilized Galerkin spatial discretization (proposed in Section 2) for the stationary advection problem. 6.1. Stationary problem: Test of convergence. Let us consider the stationary pure transport (↵ = 0) problem corresponding to (41) in the unit square ⌦ = [0, 1] 2 . We perform a set of numerical simulations on unstructured meshes {T h } h obtained by uniform refinement of an initial partition T 0 which resolves the jump discontinuity in the velocity field = ( 1 , 2 ) > . The velocity is assumed to be piecewise polynomial with respect to the open subdomain partition ⌦ 1 = (0, 0.5) ⇥ (0, 1) and ⌦ 2 = (0.5, 1) ⇥ (0, 1). Namely,

1 (x) = ⇢ 1 x 2 ⌦ 1 3 x 2 ⌦ 2 2 (x, y) = 2y + 1 in ⌦.
The data f and g in ( 41) are chosen such that the strong solution of the BVP is given by the discontinuous vector field u with components

u 1 (x) = ⇢ 3 x 2 ⌦ 1 1 x 2 ⌦ 2 u 2 (x, y) = (1 x 2 )(1 y 2 ) in ⌦.
Note that the exact solution is tangentially continuous and such that the forcing term f belongs to L 2 (⌦). We perform a numerical discretization based on:

(i) ⇤ 1 h (T h ) = P r+1 ⇤ 1 (T h ), rotated Raviart-Thomas elements (Figure 4); (ii) ⇤ 1 h (T h ) = P r ⇤ 1 (T h ), rotated Brezzi-Douglas-Marini (BDM) elements [11] (Figure 5); (iii) ⇤ 1 h (T h ) = P d r ⇤ 1 (T h
), piecewise polynomial discontinuous 1-forms (Figure 6); Figure 4, Figure 5 and Figure 6 show the behavior of the L 2 -error as the mesh is refined for the non-stabilized and stabilized Galerkin spatial scheme introduced in Section 2. As with Lipschitz continuous velocity fields, convergence rate r + 1 of the L 2 -error is attained by the stabilized scheme with edge elements of polynomial degree r. For lowest order conforming elements, the rate deteriorates by a factor of 1 when the non-stabilized scheme is applied, whereas higher order polynomial discretization yields numerical solutions which su↵er of large oscillations.

Note that a discretization based on the full polynomial space (case (ii), Figure 5) the error behaves as in the case of Lipschitz continuous velocity fields when the polynomial degree r is odd. Even polynomial degrees lead to a deteriorated convergence rate of the error in the L 2 -norm.

(T h ) = P d r ⇤ 1 (T h ).
6.2. Stationary problem: Velocity with non-resolved discontinuities. The derivation of the method, as of Section 2, relies on the assumption that the mesh resolves the possible discontinuities of the velocity field. In the following experiment we investigate an example of normal jump discontinuity in the velocity field not resolved by the mesh or any of its refinements and observe that the instabilities arising downstream of the discontinuity irremediably compromise the accuracy of the numerical solution and wreck the performance of the method. The failure of the numerical scheme in this test case may be ascribed to the fact that, since across the mesh facets the jump of the velocity vanishes, the scheme itself does not capture the discontinuity of the velocity and hence of the solution. Jump discontinuities are only taken into account through numerical quadrature.

In more details, the pure magnetic transport problem is solved in the domain ⌦ = [0, 1] 2 with a tensor product mesh and velocity field = ( 1 , 2 ) defined component-wise as

1 (x, y) = ⇢ 1 x < y 3 x > y 2 ⌘ 1 in ⌦.
The data f and g are such that the strong solution of the stationary problem corresponding to (41) is given by u = (u 1 , u 2 ) with

u 1 (x, y) = ⇢ 3 sin(⇡x) x < y sin(⇡x) x > y u 2 ⌘ (1 x 2 )(1 y 2 ) in ⌦,
as shown in Figure 7 (first column). The numerical discretization is performed with lowest order edge elements i.e. ⇤ 1 h (T h ) = Q 1 ⇤ 1 (T h ) and upwind stabilization. On the basis of Figure 7, it can be inferred that the numerical solution obtained with the upwind stabilized scheme is not a↵ected by spurious oscillations but fails to reproduce the exact solution. An error analysis provides evidence of large errors along the discontinuity and no convergence is achieved. Q 1 ⇤ 1 (T h ) and the second order Heun method as time integrator. Let S be the characteristic function on the subset S ⇢ ⌦, we consider di↵erent velocity fields whose discontinuities are resolved by the domain partition ⌦ = ⌦ 1 [ ⌦ 2 with ⌦ 1 = (0, 1) ⇥ (0, 2) and ⌦ 2 = (1, 2) ⇥ (0, 2). Numerical simulations have been conducted with the following velocity fields and initial conditions:

(i) = (0, ⌦1 ⌦2 ) > (see Figure 8 (a)) and initial condition u 0 (x, y) = (x(2 x)y(2 y), sin(2⇡x))

> : u 0 is in H(curl, ⌦), and its component in the direction of the velocity field vanishes along the discontinuity; (ii) = (0, ⌦1 ⌦2 ) > (see Figure 8 (a)) and initial condition u 0 (x, y) = (sin(2⇡x), 4)

> : The initial condition is curl-free, while its contraction with the velocity field, namely the vector component in the velocity direction is not in H > : u 0 is in H(curl, ⌦), and its component in the direction of the velocity field is not in H 1 (⌦). For the case (i), we run the simulation for an entire period, namely until T = 2 and compare the solution with the initial condition. Figure 9 shows that the initial datum is transported in the two di↵erent domains and the L 2 -error computed at the final time reaches the expected first order convergence (see Figure 10).

1 (⌦); (iii) = ( ⌦1 ⌦2 ,
)) with ⇤ 1 h (T h ) = Q 1 ⇤ 1 (T h 15
) and Heun timestepping (⌧ = 0.1h), for shear velocity and initial condition as in (i). -error

Q 1 ⇤ 1 -Heun O(h)
⇤ 1 h (T h ) = Q 1 ⇤ 1 (T h
) in space and Heun timestepping. Shear velocity and initial condition as in (i).

In case (iv), even if the initial condition is smooth, the magnetic advection with normally discontinuous collisional velocity yields the formation of a shock along the discontinuity (Figure 11 (b)) until complete blow up. An analogous behavior of the numerical magnetic potential obtained from the stabilized scheme can be reported in the case (ii), where instantaneous blow up of the solution along the discontinuity is observed. However, we expect a blow-up of the solution in these situations: the observed behavior of the numerical solution is not engendered by instabilities produced within the numerical scheme, but accurately reflects "physical reality".

A similar conclusion can be drawn for the case (iii) in Figure 12. 15)) with ⇤ 1 h (T h ) = Q 1 ⇤ 1 (T h ) and Heun method (⌧ = 0.1h), for collisional velocity and smooth initial condition as in (iv). First component of the magnetic potential obtained from the stabilized scheme (c f as in [START_REF] Crippa | Regularity and compactness for the DiPerna-Lions flow[END_REF]) with ⇤ 1 h (T h ) = Q 1 ⇤ 1 (T h ) and Heun method (⌧ = 0.1h), for collisional velocity and smooth initial condition as in (iii).

6.4. Orszag-Tang benchmark. As last numerical experiment, we rely on a widely used two dimensional benchmark problem, the so-called Orszag-Tang vortex system [START_REF] Orszag | Small-scale structure of two-dimensional magnetohydrodynamic turbulence[END_REF], which describes the transition to supersonic turbulence in the MHD equations. The development of shock waves and the complex interaction between various shocks with di↵erent speed which characterized the solution, makes the Orszag-Tang benchmark a challenging test for numerical methods.

Let us consider the two dimensional magnetic advection problem in ⌦ = [0, 2] 2 with periodic boundary conditions at the boundary @⌦. The time interval is I = [0, 1] with uniform time step ⌧ = 5• 10 3 . The initial condition is the smooth vector field u 0 (x, y) = (sin(2⇡x), sin(⇡y))

> and the velocity field is piecewise constant with respect to the mesh. In particular, it is given at each time step as the outcome of a second order Finite Volume simulation of the full MHD system (from [START_REF] Mishra | Multi-level Monte Carlo finite volume methods for nonlinear systems of conservation laws in multi-dimensions[END_REF] and [START_REF] Fuchs | Approximate Riemann solvers and robust high-order finite volume schemes for multi-dimensional ideal MHD equations[END_REF]). Note that even if the initial velocity is smooth, complex structures such as shocks and shock interactions develop in time. Concerning the discretization, the numerical scheme has been implemented on a tensor product mesh with 200 ⇥ 200 elements and H(curl, ⌦)-conforming lowest order rotated Raviart-Thomas elements ⇤ 1 h (T h ) = Q 1 ⇤ 1 (T h ) have been used for the spatial discretization, while a second order two-stage Runge-Kutta timestepping is deployed in order to exploit the mild CFL-condition of the scheme. The stabilization parameter is as in [START_REF] Crippa | Regularity and compactness for the DiPerna-Lions flow[END_REF], i.e. the upwind direction is assumed to be given by the average of the velocity field. The rotated magnetic potential obtained using the scheme [START_REF] Mishra | Multi-level Monte Carlo finite volume methods for nonlinear systems of conservation laws in multi-dimensions[END_REF] has been compared with the magnetic induction derived from the full MHD system and a second order Finite Volume discretization (from [START_REF] Mishra | Multi-level Monte Carlo finite volume methods for nonlinear systems of conservation laws in multi-dimensions[END_REF] and [START_REF] Fuchs | Approximate Riemann solvers and robust high-order finite volume schemes for multi-dimensional ideal MHD equations[END_REF]). The numerical method we have proposed well resolves shocks, is stable and no spurious oscillations occur. As can be inferred from Figure 13, the current sheet characterizing the second component of the magnetic induction (black box in Figure 13 (d)) is not captured by the H(curl, ⌦)-conforming scheme due to the low order polynomial space discretization which is highly di↵usive. Hence, if u = (! 1 , . . . , ! I , . . . ) > is a vector proxy of ! 2 ⇤ k (⌦) then P n i=1 i @ xi u + Cu with C JI := J,I ( P n i=1 @ xi i ) + e k I ^?L e k J is a vector proxy of L !.

Lemma 3 . 1 .

 31 There exists a constant C S depending on the stabilization coe cients |c f | 1/2

2 .

 2 We refer to[START_REF] Heumann | Stabilized Galerkin for Transient Advection of Di↵erential Forms[END_REF] Lemma 3.4] for a complete proof. As a consequence of the estimates shown in Theorem 3.1, Theorem 3.2 and Theorem 3.3, we present a convergence result for the stationary advection problem with Lipschitz continuous velocity fields. Analogous estimates were proposed in [23, Theorem 4.1.8] for non-conforming di↵erential forms in R n . The present result extends to H⇤ k (⌦)-conforming discrete di↵erential forms in R n the a priori convergence estimates derived in [12, Section 5] for 0-forms in R 3 and in [23, Theorem 4.1.13 and 4.1.14] for 1-and 2-forms in R 3 . Note that the numerical experiments presented in Section 6 for non-Lipschitz velocities indicate that the following convergence result might hold in a more general setting. Theorem 3.4. Let ↵ 2 L 1 (⌦) and 2 W 1,1 (⌦) in (5) satisfy the monotonicity condition (6). Furthermore, let the stabilization parameters c f fulfill the nonnegativity Assumption 2. Then

Figure 1 .

 1 Figure 1. Stabilized Galerkin discretization with upwind stabilization, polynomial di↵erential 1-forms ⇤ 1 h (T h ) and explicit Euler timestepping.

Figure 2 .

 2 Figure 2. Stabilized Galerkin discretization with upwind stabilization, polynomial di↵erential 1-forms ⇤ 1 h (T h ) and second order Heun timestepping.

  r ) (a) Non-stabilized scheme c f = 0.

  r+1 ) (b) Stabilized scheme, c f as in[START_REF] Crippa | Regularity and compactness for the DiPerna-Lions flow[END_REF].

Figure 4 .

 4 Figure 4. H(curl, ⌦)-conforming finite elements of the first family,⇤ 1 h (T h ) = P r+1 ⇤ 1 (T h ).

  r ) (a) Non-stabilized scheme c f = 0.

  r+1 ) (b) Stabilized scheme, c f as in[START_REF] Crippa | Regularity and compactness for the DiPerna-Lions flow[END_REF].

Figure 5 .

 5 Figure 5. H(curl, ⌦)-conforming finite elements of the second family, ⇤ 1 h (T h ) = Pr⇤ 1 (T h ).

  Non-stabilized scheme c f = 0.

  Stabilized scheme, c f as in[START_REF] Crippa | Regularity and compactness for the DiPerna-Lions flow[END_REF].

Figure 6 .

 6 Figure 6. Non-conforming DG finite elements ⇤ 1 h (T h ) = P d r ⇤ 1 (T h ).

6. 3 .Figure 7 .

 37 Figure 7. Vector components of the exact and numerical solutions for jumps of the velocity not resolved by the mesh.

  0) > (see Figure 8 (b)) and initial condition u 0 (x, y) = (sin(2⇡x), 4) > : The initial condition is curl-free, and its component in the direction of the velocity field vanishes along the discontinuity; (iv) = ( ⌦1 ⌦2 , 0) > (see Figure 8 (b)) and initial condition u 0 (x, y) = (x(2 x)y(2 y), sin(2⇡x))

Figure 8 .t = 2 Figure 9 .

 829 Figure 8. Sketch of shear velocity (a) and collisional velocity (b).

Figure 10 .

 10 Figure 10. L 2 -error at time T =2 of the fully discrete problem:⇤ 1 h (T h ) = Q 1 ⇤ 1 (T h) in space and Heun timestepping. Shear velocity and initial condition as in (i).

t = 1 Figure 11 .

 111 Figure 11. First component of the numerical solution obtained from the stabilized scheme (c f as in (15)) with⇤ 1 h (T h ) = Q 1 ⇤ 1 (T h) and Heun method (⌧ = 0.1h), for collisional velocity and smooth initial condition as in (iv).

t = 2 Figure 12 .

 212 Figure 12. First component of the magnetic potential obtained from the stabilized scheme (c f as in[START_REF] Crippa | Regularity and compactness for the DiPerna-Lions flow[END_REF]) with ⇤ 1 h (T h ) = Q 1 ⇤ 1 (T h ) and Heun method (⌧ = 0.1h), for collisional velocity and smooth initial condition as in (iii).

Figure 13 .

 13 Figure 13. Orszag-Tang benchmark: Comparison at time T = 1 between the magnetic potential obtained with the stabilized scheme (13)-(14), lowest order edge elements and second order Heun timestepping (first column (a)-(c)) and the corresponding components of the magnetic field obtained from the MHD simulation with a finite volume scheme in[START_REF] Fuchs | Approximate Riemann solvers and robust high-order finite volume schemes for multi-dimensional ideal MHD equations[END_REF] (second column (b)-(d)). The current sheet in Figure (d) is framed in the black box.

  Appendix A. A coordinate based representation of Lie derivatives. Let {e k I } ( n k ) I=1 be the orthonormal basis of alternating k-linear forms, e.g. ! = P ( n k ) I=1 ! I e k I TRANSIENT ADVECTION OF DIFFERENTIAL FORMS 29 with 0-forms ! I for arbitrary ! 2 ⇤ k (⌦) and e k I ^?e k J = I,J . Then the projection of L ! onto e k J yields L ! ^?e k J = L (! ^?e k J ) + ! ^?L e k J = di (! ^?e k J ) + ! ^?L e k J

Table 1 .

 1 H⇤ k (⌦) : Sobolev (Hilbert) space of k-forms, [5, Section 2.2] Notations from exterior calculus; for details see [25, Section 2] or [27, Sections 2.1 and 2.2] or compendia on di↵erential geometry.

38, 1.2.2 e)] : adjoint of the exterior derivative ? = ( 1) k d?, [38, 1.2.2 f)] i : contraction ⇤ k (⌦) ! ⇤ k 1 (⌦) with vector field [38, 1.2.2 d)] j ! : adjoint of the contraction i , ?j = ( 1) k i ?, [25, Definition 2.2 (10)] L : Lie derivative ⇤ k (⌦) ! ⇤ k (⌦) associated with vector field L : adjoint of the Lie derivative operator L tr : trace operator ⇤ k (⌦) ! ⇤ k (@⌦) [38, p. 26] ^: wedge or exterior product ⇤ k (⌦) ⇥ ⇤ `(⌦) ! ⇤ k+`( ⌦) [38, 1.2.2 a)] ? : Euclidean Hodge operator ⇤ k (⌦) ! ⇤ n k (⌦) [38, 1.2.2 c)]

Table 2 .

 2 Exterior calculus notations and corresponding expressions for vector proxies u/u of k-form !. For details see[25, Table 2], [27, Section 2.2], [5, Table 2.2].

	Exterior calculus	Proxy calculus
	^product

The term "vector proxy" was coined in[START_REF] Bossavit | On the geometry of electromagnetism. (2): 'Geometrical objects[END_REF]

partially supported by SNSF grant 146355.