High order approximation of Hodge Laplace problems with local coderivatives on cubical meshes

In mixed finite element approximations of Hodge Laplace problems associated with the de Rham complex, the exterior derivative operators are computed exactly, so the spatial locality is preserved. However, the numerical approximations of the associated coderivatives are nonlocal and it can be regarded as an undesired effect of standard mixed methods. For numerical methods with local coderivatives a perturbation of low order mixed methods in the sense of variational crimes has been developed for simplicial and cubical meshes. In this paper we extend the low order method to all high orders on cubical meshes using a new family of finite element differential forms on cubical meshes. The key theoretical contribution is a generalization of the linear degree, in the construction of the serendipity family of differential forms, and the generalization is essential in the unisolvency proof of the new family of finite element differential forms.


Introduction
In this paper we consider finite element methods for the Hodge Laplace problems of the de Rham complex where both the approximation of the exterior derivative and the associated coderivative are spatially local operators. The locality of the coderivative operator is not fulfilled in standard mixed methods for these problems (cf. [6,7]). In fact, pursuing numerical methods with local coderivative is related to the development of various numerical methods for the Darcy flow problems.
To discuss this local coderivative property in a more familiar context, let us consider the mixed form of Darcy flow problems: Find (σ, u) ∈ H(div, Ω) × L 2 (Ω) such that ∀τ ∈ H(div, Ω), div σ, v = f, v , ∀v ∈ L 2 (Ω), (1.1) where the unknown functions σ and u are vector and scalar fields defined on a bounded domain Ω in R n , and ∂Ω is its boundary. The coefficient K is symmetric, matrix-valued, spatially varying, and uniformly positive definite. Note that u is the scalar field of the pressure and σ is the fluid velocity given by the Darcy law σ = −K grad u. The standard mixed finite element method for this problem is: where Σ h ⊂ H(div, Ω) and V h ⊂ L 2 (Ω) are finite element spaces, and σ h is an approximation of σ = −K grad u. Here the notation ·, · is used to denote the L 2 inner product for both scalar fields and vector fields defined on Ω.
The mixed method (1.2) has a local mass conservation property, and its stability conditions and error estimates are well-studied (cf. [10]). However, the mixed method (1.2) does not preserve the local property of the map u → σ = −K grad u in the continuous problem. In other words, the map u h → σ h , defined by the first equation of (1.2), is not local because the inverse of the so-called "mass matrix" derived from the L 2 inner product K −1 τ, τ ′ on Σ h is nonlocal. Since constitutive laws are spatially local relations of quantities in many physical models, construction of local numerical constitutive laws is one of key issues in the development of numerical methods following physical derivation of constitutive laws such as the finite volume methods and the multi-point flux approximations.
Here we give a brief overview on previous studies of numerical methods with local coderivatives mainly for the Darcy flow problems but we have to admit that this overview and the list of literature here are by no means complete.
An early work for the locality property by perturbing the mixed method (1.2) was done in [8] on triangular meshes with the lowest order Raviart-Thomas space. The approach leads to a two-point flux method, which approximate flux across the interface of two cells by two point values of pressure field in the two cells, but the two-point flux method is not consistent in general for anisotropic K, cf. [1,2]. To circumvent this defect, various multi-point flux approximation schemes were derived (cf. [1]) but the stability and error estimates of these schemes are usually nontrivial and are restricted to low order cases. It seems that the most useful approach for the stability and error estimates for these numerical schemes is to utilize connections between the schemes and perturbed mixed finite element methods, cf. [9,16,21,22,25,27]. An alternative approach to perturbed mixed finite element methods for the local coderivative property was proposed in [11,27] independently for simplicial and quadrilateral meshes. The key to achieve local coderivatives in this approach is a mass-lumping for vector-valued finite elements. Further extensions to hexahedral grids are studied in [20,26]. Extensions to all high order methods are studied in [3] with the development of a new family of H(div) finite elements on quadrilateral and hexahedral meshes, which is inspired by the new family of low order finite elements in [23]. For the Maxwell equations, Cohen and Monk studied perturbed mixed methods based on anisotropic mass-lumping of vector-valued finite elements but the methods may not be consistent when the material coefficients are not isotropic (cf. [14]). For the Hodge Laplace problems the discrete exterior calculus, proposed in [15,19], cf. also [18], has a natural local coderivative property by construction. However, a satisfactory convergence theory seems to be limited (see [24] for 0-form case).
The purpose of this paper is to extend the results in [3,23] to the Hodge Laplace problems on cubical meshes. More precisely, we will construct high order perturbed mixed methods for the Hodge Laplace problems on cubical meshes which have the local numerical coderivatives. Since the Hodge Laplace problems are models of other specific problems, the numerical method in the present paper can be used to develop numerical methods with the locality property for other problems. For example, the methods with the newly developed finite elements have potential applications to high order mass-lumping for the time-dependent Maxwell equations. However, we will restrict our discussion only to steady problems in the paper because developing the methods for steady problems is already quite involved.
The paper is organized as follows. In the next section we will present a brief review of exterior calculus, the de Rham complex with its discretizations, and the abstract analysis results in [23] for the analysis framework to be used in later sections. In Section 3 we develop a new family of finite element differential forms on cubical meshes, sayQ r Λ k , by defining the shape functions and the degrees of freedom, and proving the unisolvency. We also prove some properties of the new elements for construction of numerical methods with the local coderivative property. In Section 4, we construct numerical methods with local coderivatives usingQ r Λ k and the framework in Section 2. Finally, we summarize our results with some concluding remarks in Section 5.

Preliminaries
Here we review the language of the finite element exterior calculus [6,7] and also introduce new concepts of polynomial differential forms. We assume that Ω ⊂ R n is a bounded domain with a polyhedral boundary. We will consider finite element approximations of a differential equation which has a differential form defined on Ω as an unknown. Let Alt k (R n ) be the space of alternating k-linear maps on R n . For 1 ≤ k ≤ n let Σ k be the set of increasing injective maps from {1, ..., k} to {1, ..., n}. Then we can define an inner product on Alt k (R n ) by where σ i denotes σ(i) for 1 ≤ i ≤ k and {e 1 , . . . , e n } is any orthonormal basis of R n . The differential k-forms on Ω are maps defined on Ω with values in Alt k (R n ). If u is a differential k-form and t 1 , . . . , t k are vectors in R n , then u x (t 1 , . . . , t k ) denotes the value of u applied to the vectors t 1 , . . . , t k at the point x ∈ Ω. The differential form u is an element of the space L 2 Λ k (Ω) if and only if the map is in L 2 (Ω) for all tuples t 1 , . . . , t k . In fact, L 2 Λ k (Ω) is a Hilbert space with inner product given by The exterior derivative of a k-form u is a (k + 1)-form du given by wheret j implies that t j is not included, and ∂ tj denotes the directional derivative. The Hilbert space HΛ k (Ω) is the corresponding space of k-forms u on Ω, which is in L 2 Λ k (Ω), and where its exterior derivative, du = d k u, is also in L 2 Λ k+1 (Ω). The L 2 version of the de Rham complex then takes the form In the setting of k-forms, the Hodge Laplace problem takes the form where d = d k is the exterior derivative mapping k-forms to (k + 1)-forms, and the coderivative d * = d * k can be seen as the formal adjoint of d k−1 . Hence, the Hodge Laplace operator L above is more precisely expressed as A typical model problem studied in [6,7] is of the form (2.2) and with appropriate boundary conditions. The mixed finite element methods are derived from a weak formulation, where σ = d * u is introduced as an additional variable. It is of the form: Here ·, · denotes the inner products of all the spaces of the form L 2 Λ j (Ω) which appears in the formulation, i.e., j = k−1, k, k+1. We refer to Sections 2 and 7 of [6] for more details. We note that only the exterior derivate d is used explicitly in the weak formulation above, while the relation σ = d * k u is formulated weakly in the first equation. The formulation also contains the proper natural boundary conditions. The problem (2.3) with k = n − 1 corresponds to a weak formulation of the elliptic equation (2.2) in the case when the coefficient K is the identity matrix. The weak formulations (2.3) can be modified for variable by changing the L 2 inner products to the inner product with the variable coefficient, see [6,Section 7.3]. Throughout the discussion below we will restrict the discussion to the constant coefficient case but the extension of the discussion to problems with variable coefficients which are piecewise constants with respect to the mesh, is straightforward. We refer to [23,Section 6] for details.
If the domain Ω is not homologically trivial, then there may exist nontrivial harmonic forms, i.e., nontrivial elements of the space H k (Ω) = {v ∈ HΛ k (Ω) : dv = 0 and v, dτ = 0 for all τ ∈ HΛ k−1 (Ω)}, and the solutions of the system (2.3) may not be unique. To obtain a system with a unique solution, an extra condition requiring orthogonality with respect to the harmonic forms: The basic construction in finite element exterior calculus is a corresponding subcomplex of (2.1), where the spaces V k h are finite dimensional subspaces of HΛ k (Ω). In particular, the discrete spaces should have the property that d( The finite element methods studied in [6,7] are based on the weak formulation (2.3). These methods are obtained by simply replacing the Sobolev spaces HΛ k−1 (Ω) and HΛ k (Ω) by the finite element spaces V k−1 h and V k h . More precisely, we are searching for a triple h , approximating the harmonic forms, is given by Stability and error estimates for the numerical solution of (2.5) are discussed in [7,Theorem 3.9] with an error estimate of the form with (σ, u, p) X := σ 2 + dσ 2 + u 2 + du 2 + p 2 1 2 and E h (u) comes from the nonconformity of the space of discrete harmonic forms, H k h , to the space of continuous harmonic forms H k . E h (u) vanishes if there is no nontrivial harmonic forms, and will usually be of higher order than the other terms on the right-hand side of (2.6). We refer to [6, Section 7.6] and [7, Section 3.4] for more details.
In [23], abstract conditions are proposed to develop numerical methods satisfying such properties and the lowest order methods were developed in simplicial and cubical meshes. In particular, the cubical mesh case required construction of new finite element differential forms which were called S + 1 Λ k spaces. In [3], the S + 1 Λ k family is extended to all higher orders for the Darcy flow problems in the two and three dimensions on cubical meshes, and as a consequence, higher order numerical methods satisfying the local coderivative property in the two and three dimensional Darcy flow problems on weakly distorted quadrilateral and hexahedral meshes are constructed. The goal of this paper is developing high order numerical methods for the Hodge Laplace equations satisfying the local coderivative property on cubical meshes. For the stability and error analysis we use the abstract framework in [23]. The main contributions of this paper are construction of a new family of finite element differential forms, a higher order extensions of S + 1 Λ k , and the new techniques for the development of the new elements. Although the new family is a higher order extension of S + 1 Λ k , we will call itQ r Λ k family because it is closely related to the Q r Λ k space rather than the S r Λ k family. One characteristic feature ofQ r Λ k family is that the number of degrees of freedom is same as the one of Q r Λ k .
Here we summarize the abstract conditions for stability and error estimates established in [23].
and a linear map Π h : If these assumptions are satisfied, then we find a solution (σ h , u h , p h ) of the problem The following result is proved in [23].

Construction ofQ r Λ k
In this section we construct a new family of finite element differential forms Q r Λ k (T h ) on cubical meshes T h of Ω where the elements in T h are Cartesian product of intervals. If r = 1, thenQ r Λ k is the S + 1 Λ k space in [23]. For n = 2, 3 and k = n − 1,Q r Λ k spaces are the H(div) finite element spaces which were discussed in [3]. The idea of these new elements construction in [23] and [3] is enriching the shape function space of the Q − 1 Λ k and the Raviart-Thomas-Nedelec elements on cubical meshes with d-free and divergence-free shape functions in order to associate a basis of the shape function space to the degrees of freedom given by nodal point evaluation. The most technical difficulty of new element construction in [3] is the unisolvency proof. For this, the authors in [3] used some features of the enriched shape functions observed from their explicit expressions. However, it is highly nontrivial to use the same argument toQ r Λ k for general n and k because explicit forms of the enriched shape functions for general n and k are too complicated to use conventional unisolvency arguments. To circumvent this difficulty we will introduce new quantities of polynomial differential forms which allow us to extract useful features of polynomial differential forms for unisolvency proof.
In the paper the space of polynomial or polynomial differential forms without specified domain will denote the space on R n . For example, Q r Λ k is the space of polynomial coefficient differential forms such that the polynomials coefficients are in Q r (R n ). We use PΛ k to denote the space of all differential forms with polynomial coefficients.
For later discussion we introduce some additional notation for cubes in a hyperspace in R n . For given I = {i 1 , . . . , i m } ⊂ {1, . . . , n} with i 1 < i 2 < . . . < i m , consider an m-dimensional hyperspace F in R n determined by fixing all x j coordinates of j ∈ I. If f is an m-dimensional cube in the m-dimensional hyperspace F , then Σ k (f ) with k ≤ m is defined by the set of increasing injective map from {1, . . . , k} to I. For σ ∈ Σ k (f ) we will use σ to denote the range of σ, i.e., where the coefficients u σ 's are scalar functions on Ω. Furthermore, the exterior derivative du can be expressed as is defined by contraction with the vector x, i.e., (κu) x = x u x . As a consequence of the alternating property of Alt k (R n ), it therefore follows that κ • κ = 0. It also follows that where dx σi means that the term dx σi is omitted. This definition is extended to the space of differential k-form on Ω by linearity, i.e., For future reference we note that If f is an (n−1)-dimensional hyperspace of R n obtained by fixing one coordinate, for example, For a multi-index α of n nonnegative integers, x α = x α1 1 · · · x αn n . If u is in Q r Λ k , the space of polynomial k-forms with tensor product polynomials of order r, then u can be expressed as Denoting H r Λ k the space of differential k-forms with homogeneous polynomial coefficients of degree r, we also have the identity cf. [6,Section 3]. Finally, throughout this paper we setT = [−1, 1] n .
3.1. The shape function space and the degrees of freedom ofQ r Λ k . In this subsection we define the shape function space ofQ r Λ k . The shape functions of Q r Λ k will be obtained by enriching the shape functions of Q − r Λ k . There are other families of cubical finite element differential forms, cf. for example [4,12,13,17], but these spaces are not involved in the construction of our new elements.
If m is a k-form given by m = p dx σ , where σ ∈ Σ k and the coefficient polynomial p is a monomial, then we will call m form monomial. Before we define the shape functions ofQ r Λ k let us introduce new quantities of form monomials. For a poly- , resp.) conforming indices (nonconforming indices, resp.). For a form monomial m = c α x α dx σ = 0 and α i = s for a conforming (nonconforming, resp.) index i of m, we call i a conforming (nonconforming, resp.) index of degree s. We also define the conforming and nonconforming s-degrees of m = c α x α dx σ by We can define cdeg s (v) and ncdeg s (v) if these quantities are same for any form monomial of v. We remark that ncdeg 1 (m) is same as the linear degree in [4] for a form monomial m. In contrast to the linear degree, we do not define cdeg s and ncdeg s for general polynomial differential forms.
The conforming degree gives a new characterization of the shape function space (3.5) it is easy to see that We define the shape function spaceQ r Λ k as Lemma 3.1. Let m = x α dx σ ∈ PΛ k , 1 ≤ k ≤ n − 1, be given with a multi-index α and a positive integer s. Assume that m ′ and m ′′ are given as ]. Then the following identities hold: where δ i,j is the Kronecker delta. In particular, if m ∈ B r Λ k and m ′′ is a form monomial of κm with j which is a nonconforming index of degree (r + 1), then j is a conforming index of degree r in m. For the particular case, if m ∈ B r Λ k and m ′′ is a form monomial of κm which has a nonconforming index of degree (r + 1), then the nonconforming index of m ′′ must be σ j and α σj = r because α l ≤ r for 1 ≤ l ≤ n and α σj + 1 = r + 1. This completes the proof. Proof. It is easy to check the assertions by Lemma 3.1.
We now prove thatQ r Λ k is invariant under dilation and translation.
Lemma 3.4. If φ : R n → R n is a composition of dilation and translation, then Proof. Let φ(x) = Ax + b for a given invertible n × n diagonal matrix A and a vector b ∈ R n . To show φ * Q r Λ k ⊂Q r Λ k , assume that u ∈Q r Λ k is written as where we used φ * κu + = κφ * u + + b (φ * u + ) in the last equality (cf. [6, Section 3.2]). We can easily check φ * u − ∈ Q − r Λ k from the definition of φ * , and dκQ − r Λ k ⊂ Q − r Λ k from (3.9) and (3.8). From (3.6) we have Lemma 3.5. The following hold: (a) For a form monomial m = 0 in B r Λ k , κm generates at least one form monomial whose nonconforming (r + 1)-degree is 1.
Proof. For m = x α dx σ ∈ B r Λ k there is at least one conforming index σ i such that α σi = r. Then x α κ(dx σ ) has at least one form monomial such that σ i is a nonconforming index and its polynomial coefficient has x r+1 σi as a factor, so (a) is proved.
For the injectivity of dκ on B r Λ k , it suffices to show that κ is injective on B r Λ k because d is injective on the image of κ. To show κ is injective on B r Λ k , we show that form monomials with positive nonconforming (r + 1)-degree generated by κm for m ∈ B := {x α dx σ : cdeg r (x α dx σ ) > 0} are distinct. More precisely, if κm and κm for m,m ∈ B have a same form monomial (up to ±1) whose nonconforming (r + 1)-degree is 1, then m =m. To show it by contradiction, let m = x α dx σ andm = xαdxσ be two distinct elements in B and assume that κm and κm have a common form monomial with nonconforming index of degree (r + 1). From the definition of κ and the common form monomial assumption, there exist Since σ i andσ i are the only nonconforming indices of degree (r + 1) by Lemma 3.1, σ i =σ i and therefore dx σ = dxσ. Moreover, comparison of x α and xα leads to α =α, so it contradicts to x α dx σ = xαdxσ.
The following key result is a consequence of the above lemma.
Proof. By Lemma 3.5 (b), the spaces dκB r Λ k and B r Λ k have the same dimension, so it suffices to show that Q − r Λ k ∩ dκB r Λ k = {0}. Suppose that 0 = u ∈ B r Λ k and dκu ∈ Q − r Λ k . Every form monomial m of dκu satisfies cdeg r (m) ≥ 1 by (3.5) and Lemma 3.1 because ncdeg r+1 (m) = 0. However, it is a contradiction to the characterization of Q − r Λ k in (3.4), so u = 0. 3.2. Degrees of freedom and unisolvence ofQ r Λ k . In this subsection we define the degrees of freedom ofQ r Λ k and prove the unisolvency for the degrees of freedom.
For the degrees of freedom we define two polynomial spaces for σ ∈ Σ k (f ) for a cube f included in an m-dimensional hyperspace by where P r (x i ) is the space of polynomials of x i with degree less than or equal to r.

Theorem 3.7. The number of degrees of freedom given by (3.17) is
We can easily check that Therefore the number of degrees of freedom given by (3.17) is so the proof is complete.
The following result will be useful to reduce unisolvency proof.
Theorem 3.8 (trace property). Let f be a hyperspace in R n determined by fixing one coordinate x i . Then Proof. Without loss of generality we assume that f = {x ∈ R n : x n = c} for some constant c. It is easy to check by definition that tr f Q − r Λ k ⊂ Q − r Λ k (f ), so it is enough to show that tr f dκB r Λ k ⊂Q r Λ k (f ). We will show tr f dκ(x α dx σ ) ∈ Q r Λ k (f ) for all x α dx σ ∈ B r Λ k below. By (3.2) and the commutativity of tr f and d, , then x f (x α dx σ ) = 0. Since tr f (x α dx σ ) = x α | xn=c dx σ ∈ Q r Λ k (f ), the conclusion follows from (3.14).
Before we start the unisolvence proof we need a lemma and auxiliary definitions.
Proof. We first prove it under the assumption i >ĩ. Let The proof is complete.
By Lemma 3.1 dκ maps D r,l Λ k into itself. Considering the decomposition of B r Λ k Let us recall the degrees of freedom of Q − r Λ k (T ) with vanishing trace in [5]. If u ∈ Q − r Λ k (T ) and tr f u = 0 for all f ∈ ∆ n−1 (T ), then implies that u = 0.
We are now ready to prove unisolvency ofQ r Λ k (T ), r ≥ 1, with the degrees of freedom (3.17).
Proposition 3.10 (unisolvence with vanishing trace assumption). Suppose that u ∈Q r Λ k (T ), r ≥ 1, and tr f u = 0 for all f ∈ ∆ n−1 (T ). If Proof. If k = 0 or k = n, thenQ r Λ k (T ) = Q r Λ k (T ) and (3.20) is a standard set of degrees of freedom for the shape functions with vanishing trace, so there is nothing to prove.
Assume that 0 < k < n, and let u = σ∈Σ k (T ) u σ dx σ ∈Q r Λ k (T ) be a shape function with vanishing trace. From the vanishing trace assumption tr f u = 0 for all f ∈ ∆ n−1 (T ), u σ vanishes on all faces f ∈ ∆ n−1 (T ) determined by l ) for all σ ∈ Σ k (T ). In the degree properties (3.8), (3.9), (3.10), the coefficients of all form monomials inQ r Λ k (T ) have at most one variable of degree (r + 1). From this observation and (3.20), u has a form where L w s (t) is the monic Legendre polynomial of degree s on [−1, 1] with weight (1 − t 2 ), and p σ,i is a polynomial in (Q r,σ ⊗ Q r−2,σ * )(T ) which is independent of x i . Recall the decomposition (3.18) and let Let l 0 be the largest index l such that u l,s = 0 for some s, and let s 0 be the largest index s such that u l0,s = 0. In the proof below, we will show u l0,s0 = 0. By induction, this implies that u = u 0 ∈ Q − r Λ k (T ), and therefore u = 0 by (3.20), (3.19), and the inclusion We now start the proof of u l0,s0 = 0. In the form (3.21), if p 0 is a monomial of p σ,i for fixed σ, then cdeg r (p 0 dx σ ) < l 0 because any form monomial m containing the factors x r+1 i and p 0 in its polynomial coefficient, which comes from the expansion of the coefficient of dx σ in (3.21), satisfies both of ncdeg r+1 (m) = 1 and ncdeg r+1 (m) + cdeg r (m) ≤ l 0 . Note that it also implies that u l0,s0 does not have any form monomials with conforming r-degree l 0 .
We now show u l0,s0 = 0 for l 0 = 1. If l 0 = 1, then q σ,i ∈ (Q r−1,σ ⊗ Q r−2,σ * )(T ) because cdeg r (q σ,i dx σ ) < l 0 . We claim that q σ,i = 0 for all σ and i. To prove it we show that the form monomials in du l0,s0 with conforming r-degree 1, are all distinct. Note that such form monomials are from . For fixed σ the form monomials in this formula are all distinct for different i's. Further, if we assume that there isσ ∈ Σ k (T ) withσ = σ andĩ ∈ [[σ * ]] such that the form monomials in are not linearly independent with the ones in (3.23), then it leads to a contradiction because i andĩ are the only conforming indices of degree r in these form monomials, and therefore i =ĩ, σ =σ. Thus, all form monomials in du l0,s0 with conforming r-degree 1 are distinct and have the form (3.23). From du l0,s0 = 0, q σ,i = 0 for all σ ∈ Σ k (T ) and i ∈ [[σ * ]], and therefore u l0,s0 = 0 by (3.22). If l 0 > 1, then we consider the expression of q σ,i as a sum of homogeneous polynomials in which q σ,i,τ ∈ (Q r−1,σ ⊗ Q r−2,σ * )(T ) is independent of the variables x i and x j 's for j ∈ τ . Here q σ,i,τ can be vanishing for some τ ⊂ [[σ]]. In the discussion below, we will make a system of equations with all possibly involving q σ,i,τ 's as its unknown. The equations are obtained from κu l0,s0 and du l0,s0 , and the right-hand sides are zero. We then show that the system is well-posed, so all q σ,i,τ 's vanish.
To derive equations from κu l0,s0 , note that κu l0,s0 ∈ κdκB r Λ k (T ) = κB r Λ k (T ), so the nonconforming (r + 1)-degree of form monomials in κu l0,s0 is at most 1 by Lemma 3.1. In the formal expression of κu l0,s0 using (3.22) and (3.24), for fixed σ, the form monomials which have nonconforming (r + 1)-degree 2 are obtained by ] with |τ |= l 0 − 1, and j ∈ τ . We claim that all of these terms are linearly independent for fixed σ ∈ Σ k (T ). To see it, assume that the form (3.25) with (i, τ, j) and (i ′ , τ ′ , j ′ ) are linearly dependent for Comparison of the nonconforming indices with degree (r + 1) gives i = i ′ , j = j ′ or i = j ′ , j = i ′ . If i = i ′ and j = j ′ , then τ = τ ′ from the comparison of conforming indices of degree r. If i = j ′ and j = i ′ , then dx σ−j = dx σ−j ′ , so they cannot be linearly dependent. Therefore the terms of the form (3.25) are all distinct for fixed σ.
Proof. If l = k, then tr f u ∈Q r Λ k (f ) = Q r Λ k (f ) for f ∈ ∆ k (T ) by Theorem 3.8, and (3.30) gives a standard set of degrees of freedom for Q r Λ k (f ), so tr f u = 0.
Suppose that tr f u = 0 for all f ∈ ∆ l (T ) for some l ≥ k. For any givenf ∈ ∆ l+1 (T ), trf u ∈Q r Λ k (f ) and tr f trf u = tr f u = 0 for all f ∈ ∆ l (f ) by the assumption. By Proposition 3.10 and the degrees of freedom (3.30), trf u = 0. By induction, we can show that tr f u = 0 for all f ∈ ∆ n−1 (T ), so u = 0 by Proposition 3.10.

3.3.
Nodal tensor product degrees of freedom. In this subsection we show thatQ r Λ k (T ) has a set of degrees of freedom given by evaluating nodal values at a set of points inT . This alternative set of degrees of freedom will be used to define numerical methods with local coderivatives in the next section.
The Gauss-Lobatto quadrature rule with r + 1(r ≥ 1) points uses r − 1 interior points and two end points of I = [−1, 1] with positive weights and gives exact integration of polynomials of order 2r − 1. Let {ξ j } r j=0 be the quadrature points of the Gauss-Lobatto quadrature on I and {λ j } r j=0 with λ j > 0 be the weights at the points. We can define a quadrature rule onT by taking tensor products of the Gauss-Lobatto quadrature nodes and weights, and we discuss the formal expressions for the tensor product quadrature rules below.
For σ ∈ Σ k (T ) let N σ,r be the set of points in 1≤i≤k I(x σi ) defined by N σ * ,r is defined similarly as N σ * ,r = {(x σ * 1 , . . . , x σ * n−k ) = (ξ j1 , . . . , ξ j n−k ) : 0 ≤ j l ≤ r, 1 ≤ l ≤ n − k}, the set of points in the (n − k)-dimensional cube 1≤i≤n−k I(x σ * i ). Therefore the tensor product N r := N σ,r ⊗ N σ * ,r is the set given by the tensor product of Gauss-Lobatto quadrature points onT . Letting N be the set of nonnegative integers, we will use ξ i σ and ξ j σ * with multi-indices i ∈ N k and j ∈ N n−k to denote the points in N σ,r and N σ * ,r , respectively. We also use λ j σ and λ j σ * to denote the corresponding tensor-product weights of the Gauss-Lobatto quadrature.
For a continuous function v onT , we define Similarly, for polynomial differential forms u = u σ dx σ ∈ PΛ k (T ) we define R r (u) as the element in R M ⊗ Λ k with M = n k (r + 1) n , which consists of R r (u σ ) ⊗ dx σ for σ ∈ Σ k (T ).
Suppose also that nonconforming (r + 1)-degree of every form monomial in v σ dx σ is at most 1 for all σ ∈ Σ k (T ).
Proof. If we rewriteṽ σ with the weighted Legendre polynomial L w j 's, then the assumption on nonconforming (r + 1)-degree leads us to havẽ for all i ∈ [[σ * ]] and ξ i σ ∈ N r,σ . To see it, note that the quadrature along x i coordinate can be replaced by integration with x i variable on [−1, 1] because the degree of x i variable is r If we consider the expression ofṽ σ,0 in (3.32), a completely analogous argument using orthogonality of Legendre polynomials gives . Note that this result is obtained without using the assumption R r (v σ ) = 0. Since R r (v σ ) = 0, the above quantity vanishes. By the definition of E r,ξ i σ , we have either b σ * (ξ j σ * )ψ(ξ j σ * ) = 0 for all ξ j σ * ∈ N r,σ * or q d ′ (ξ i σ ) = 0. However, the first case implies that b σ * ψ = 0 because the nodal value evaluations at the points in N r,σ * is already a set of degrees of freedom for Q r,σ * (T ), and it leads to a contradiction to ψ = 0. Therefore q d vanishes at ξ i σ , and we can show that q d vanishes at any point in N r,σ with the same argument. Recall that q d ′ ∈ Q r,σ (T ), so these vanishing conditions of q d ′ implies that q d ′ = 0. Finally, this holds for any d ′ , and thereforeṽ σ,0 = 0.
To show v = 0, we notice that v σ with (3.31) is exactly the form of u σ in (3.21), so the same argument in the unisolvency proof can be used to show v = 0. Proof. Note that tr f v ∈Q r Λ k (f ) = Q r Λ k (f ) for f ∈ ∆ k (T ). Since the restriction of R r on f becomes a set of quadrature degrees of freedom of Q r Λ k (f ), tr f v = 0 for all f ∈ ∆ k (T ) holds. For tr g v with g ∈ ∆ k+1 (T ), all traces of tr g v on kdimensional subcubes are vanishing, so the assumption of Lemma 3.12 is satisfied for tr g v. Applying Lemma 3.12 with the restriction of R r on g, one can conclude that tr g v = 0 for any g ∈ ∆ k+1 (T ). We can continue this argument for any g ∈ ∆ l (T ), l ≥ k + 1, by induction, so the conclusion follows.

Numerical methods with local coderivatives
We construct numerical methods with local coderivatives usingQ r Λ k (T h ). For this we need a modified bilinear form ·, · h in (A), and the auxiliary spacesṼ is the space of discontinuous polynomial differential forms. We will show that the finite elementsQ r Λ k (T h ) and Q − r Λ k (T h ) satisfy the conditions (A), (B). In addition to the conditions in (A) and (B), for local coderivatives, we also need to show that ·, · h onQ r Λ k (T h ) can give a block diagonal matrix with appropriate choice of global DOFs.
We have shown thatQ r Λ k (T ) has a set of DOFs determined by evaluations at nodal points. For T ∈ T h we can define the evaluation operator E T r for continuous functions on T with the scaled Gauss-Lobatto quadrature rules on T . For It is easy to see that ·, · h,T is an inner product onQ r Λ k (T ) and the norm defined by this inner product is equivalent to the L 2 norm with constants independent on the scaling of T . We can define the inner product ·, · h onQ r Λ k (T h ) by for u, v ∈Q r Λ k (T h ), and the norm · h defined by this inner product is equivalent to the L 2 norm with constants independent of h. Therefore (A) is satisfied.
We now verify (B), but with k-forms instead of (k − 1)-forms for notational convenience. Recall that we already definedṼ k h and W k h in (4.1), so verification of (B) consists of the following three steps: Step 1. Prove (2.7) Step 2. Define Π h satisfying the conditions in (B) Step 3. Prove (2.8)

For
Step 1, it is enough to show the identity in (2.7) for the restrictions of polynomial differential forms on any T ∈ T h . Moreover, by scaling argument, it Before we start its proof, recall that the quadrature nodes of the Gauss-Lobatto rule with (r + 1) points are the zeros of d dt L r (t) on [−1, 1]. We also note that (r + 1)(L r+1 (t) − L r−1 (t)) = (2r + 1)(tL r (t) − L r−1 (t)) = (2r + 1) i.e., the evaluation of L r+1 (t) − L r−1 (t) at the quadrature nodes of the Gauss-Lobatto rule with (r + 1) points vanishes.
As Π h , we defineΠ h : where p σ,i ∈ Q σ * ,r (T )⊗Q σ,r−1 (T ), and every term in v σ,0 written with the Legendre polynomials has a factor of L r (x j ) for some j ∈ [[σ]].
To prove the assertion by induction we need to show the following two results. First, for f ∈ ∆ k (T ) every polynomial coefficient of tr f v ∈Q r Λ k (f ) satisfies (4.7). Second, if every polynomial coefficient of tr f v satisfies (4.7) for all f ∈ ∆ l (T ) with given l ≥ k, then tr g v satisfies the same property corresponding to (4.7) for all g ∈ ∆ l+1 (T ).
We prove the first claim. Recalling tr f v ∈Q r Λ k (f ) = Q r Λ k (f ) for f ∈ ∆ k (T ) and the definition ofΠ h , the polynomial coefficient of tr f v is a polynomial in Q r (f ) which is orthogonal to all Q r−1 (f ) in the L 2 inner product. Then it has at least one factor of L r (x j ) for j ∈ [[σ]], therefore it is a form of (4.7).
For the second claim, from the trace property and the definition ofΠ h commuting with the trace operator tr, we can reduce the claim to the case l = n− 1 without loss of generality. Suppose that every polynomial coefficient of tr f v satisfies (4.7) for all f ∈ ∆ n−1 (T ).
By (4.4) and the definition of shape function space ofQ r Λ k (T ), v σ has a form where p σ,i ∈ (Q σ * ,r ⊗ Q σ,r−1 )(T ) but p σ,i is independent of x i , v σ,0 is as in the assertion, and v σ,1 is a polynomial such that every term in its expression with the Legendre polynomials, has a factor of L r (x l ) or L r−1 (x l ) for some l ∈ [[σ * ]]. We remark that v σ,0 may have polynomial terms of degree greater than r. Let us where q σ ∈ (Q σ * ,r ⊗ Q σ,r−1 )(T ). Note that q σ is a polynomial such that every term in its Legendre polynomial expression has L r (x l ) factor for some l ∈ [[σ * ]] or is in (L r+1 (x l ) − L r−1 (x l )) p σ,l | xi=1 +v σ,0 | xi=1 +q σ | xi=1 .
Proof. It suffices to show the assertion for u ∈Q r Λ k (T ) and v ∈ Q d r−1 Λ k (T ) for any T ∈ T h . We defineû ∈Q r Λ k (T ) andv ∈ Q d r−1 Λ k (T ) as the pullbacks φ * u and φ * v with the dilation map φ :T → T . From the definition of Π h and the property of pullback on differential forms, φ * (u − Π h u) =û −Π hû holds. Since u − Π h u, v h,T is a constant multiple of û −Π hû ,v h,T with the constant depending on scaling.
Finally, we check that the proposed method give local coderivatives with an argument completely analogous to the one in [23].
To see more details, let the set of nodal degrees of freedom and the basis of Q r Λ k (T h ) associated to the nodal degrees of freedom be {φ (T,z) : z ∈ N r (T ), T ∈ T h }, {ψ (T,z) : z ∈ N r (T ), T ∈ T h }, such that φ (T,z) (ψ (T ′ ,z ′ ) ) = δ T T ′ δ zz ′ with the Kronecker delta. For each z ∈ N r (T ) for some T ∈ T h and we can consider the set of all basis functions ψ (T ′ ,z) with the common point z. We denote this set of basis functions by Ψ z . Then the quadrature bilinear form ·, · h with {ψ (T,z) } gives a block diagonal matrix such that each block is associated to Ψ z for some z. The influence of the inverse of this block matrix associated to Ψ z is only on the local domain consists of {T ′ } such that ψ (T ′ ,z) ∈ Ψ z , so it is a spatially local operator and then gives a local coderivative.

Concluding remarks
In this paper we develop high order numerical methods with finite elements for the Hodge Laplace problems on cubical meshes that admit local approximations of the coderivatives. The main task of the development is construction of a new family of finite element differential forms on cubical meshes which satisfy the properties for the stability analysis framework presented in Section 2.
In our study we have only considered Hodge-Laplace problems with the identity coefficients. However, the discussion can be easily extended to a matrix-valued coefficient K −1 which is symmetric positive definite, and piecewise constant on T h . We point out that this is an important difference between our approach and the work in [14] for the Maxwell equations. More precisely, the quadrature rules in [14] are combinations of Gauss and Gauss-Lobatto rules for different components, so the methods are not robust for matrix-valued coefficients because the matrixvalued coefficients in bilinear forms mix the components with different quadrature rules. We refer the interested readers to [23,Section 6] for details of the analysis with K −1 = I.
Finally, the extension of the methods to curvilinear meshes is not easy because the conditions in (B) strongly rely on the properties of quadrature rules on cubical meshes which do not hold on distorted cubical (or called curvilinear) meshes. In [3], this extension was done for weakly distorted quadrilateral and hexahedral meshes by measuring changes of bilinear forms and related quantities under mesh distortion. However, this is possible in the cases because explicit forms of the bilinear and trilinear maps and their derivatives are available for the specific setting k = n − 1 and n = 2, 3. The extension to general n and k is completely open.