Marien-Lorenzo Hanot
email: marien-lorenzo.hanot@umontpellier.fr

AN ARBITRARY ORDER AND POINTWISE DIVERGENCE-FREE FINITE ELEMENT SCHEME FOR THE INCOMPRESSIBLE 3D NAVIER-STOKES EQUATIONS

Keywords: Finite Element, exterior calculus, incompressible Navier-Stokes, Hodge decomposition, de Rham complex, mixed element AMS subject classifications. 35Q30, 65N30, 76D07, 76M10

In this paper we discretize the incompressible Navier-Stokes equations in the framework of finite element exterior calculus. We make use of the Lamb identity to rewrite the equations into a vorticity-velocity-pressure form which fits into the de Rham complex of minimal regularity. We propose a discretization on a large class of finite elements, including arbitrary order polynomial spaces readily available in many libraries. The main advantage of this discretization is that the divergence of the fluid velocity is pointwise zero at the discrete level. This exactness ensures pressure robustness. We focus the analysis on a class of linearized equations for which we prove well-posedness and provide a priori error estimates. The results are validated with numerical simulations.

Introduction.

Recently much work has been done to design structure preserving methods, but while the construction of such methods was found early on in two dimensions, the three-dimensional case remained difficult and the introduction of the finite element exterior calculus brought a significant breakthrough. An excellent review is given by V. John et al. [START_REF] John | On the Divergence Constraint in Mixed Finite Element Methods for Incompressible Flows[END_REF]. The general idea taken from the finite element exterior calculus is to use a subcomplex of the De Rham complex. There are well-known discrete counterparts of this complex with minimal regularity, however the discretization of smoother variants is still an active topic usually leading to shape functions of high degree, for example see [START_REF] Neilan | Discrete and conforming smooth de Rham complexes in three dimensions[END_REF]. We chose to use the complex with minimal regularity, as this is often done for electromagnetism or recently for magnetohydrodynamics (see [START_REF] Hu | Helicity-conservative finite element discretization for incompressible MHD systems[END_REF]).

The main difference from usual schemes lies in the regularity of the velocity field since we only require it to be in H(div) and in the discrete adjoint of H(curl). Although the continuous space regularity is the same as the usual one, since the adjoint of (curl, H(curl)) is (curl, H 0 (curl)), and the velocity is sought in H(div) ∩ H 0 (curl) ⊂ H 1 (for a smooth enough domain, it is discussed in part 3.2 of [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF]). This does not hold (in general) in the discrete case since for V h ⊂ H(curl) and (curl * , V * h ) the adjoint of (curl, V h ) we no longer have V * h ⊂ H(curl). This has a fundamental impact both from the philosophical and practical point of view. In practice v ∈ H(div) (resp. v ∈ H(curl)) does not impose continuity of the tangential (resp. normal) components on faces. This suggests that we will not have any degree of freedom corresponding to these and lack any way to set them in a Dirichlet boundary condition. This means that the normal and the tangential part of the boundary condition must be treated in two different ways. It makes more sense in the exterior algebra and means that the fluid velocity is really sought as a 2-form (mostly defined by its flux across cell boundaries) which happens to be regular enough to also be in the domain of the exterior derivative adjoint.

Let us summarize the main idea of the algorithm. In order to preserve the free divergence constraint, we have to consider u as a 2-form, which can be discretized by face elements. Then it is not straightforward to discretize the Laplacian in the usual way ∇u, ∇v because ∇u is not a natural quantity for a 2-form. Our simple trick is to use ∇ • u = 0 to rewrite the Laplacian:

(1.1) ∆u =∇(∇ • u) -∇ × (∇ × u) = -∇ × (∇ × u).
Let Ω be a bounded domain of R 3 and T > 0, we recall the Navier-Stokes equations:

(1.2) u t + (u • ∇)u -ν∆u + ∇p = f on Ω × (0, T ),

∇ • u = 0 on Ω × (0, T )
together with some boundary and initial conditions, where u is the velocity of the fluid, p the pressure, ν the kinematic viscosity and f an external force.

Using the Lamb identity (u • ∇)u = (∇ × u) × u + 1 2 ∇(u • u) and Equation (1.1), we get the following formulation (1.3) u t + (∇ × u) × u + ν∇ × (∇ × u) + ∇P = f on Ω × (0, T ), ∇ • u = 0 on Ω × (0, T ).

where P = p + 1 2 u • u is the Bernoulli pressure. Since u is a 2-form it is not natural to take ∇ × u (and it is unadvisable for reasons detailed in Remark 3.3). Therefore, we introduce an auxiliary variable ω = ∇ × u (namely the vorticity) and work with a mixed problem. This is known as the vorticityvelocity-pressure formulation and was considered by many others (see [START_REF] Arnold | Mixed finite element approximation of the vector Laplacian with Dirichlet boundary conditions[END_REF][START_REF] Amoura | Spectral element discretization of the vorticity, velocity and pressure formulation of the Navier-Stokes problem[END_REF][START_REF] Kreeft | Mixed mimetic spectral element method for Stokes flow: A pointwise divergence-free solution[END_REF][START_REF] Dubois | Vorticity-velocity-pressure and stream function-vorticity formulations for the Stokes problem[END_REF][START_REF] Anaya | Analysis and Approximation of a Vorticity-Velocity-Pressure Formulation for the Oseen Equations[END_REF]). The finite element exterior calculus framework is very flexible because it allows us to work in abstract spaces which can be discretized easily provided that some exactness properties are fulfilled. Therefore we shall use a generic name for our spaces, here V 1 × V 2 × V 3 for the continuous spaces and V 1 h × V 2 h × V 3 h for the discrete spaces (indexed by the mesh size h). The exact requirements for these spaces require introducing some concepts and notations, hence for the sake of readability we postpone the definition to Section 2.1. Typically a valid choice is to take V 1 = H(curl, Ω), V 2 = H(div, Ω) and V 3 = L 2 (Ω). We also need another space H 3 ⊂ V 3 . This is a vector space that does not only depend on V 3 but on the couple V 2 × V 3 and which is typically of small dimension, i.e. of dimension 0 or 1 in the case of interest to us. From the point of view of exterior calculus H 3 is the set of harmonic 3-forms, for a more practical point of view H 3 can be understood as a set of Lagrange multipliers. 3 . More details are given later in Section 2.2. An example of discrete in time, mixed and linearized weak formulation is: Given

Remark 1.1. The choice of boundary conditions is encoded in the choice of V 1 × V 2 × V
f n ∈ L 2 (Ω), find (ω n , u n , p n , φ n ) ∈ V 1 × V 2 × V 3 × H 3 such that ∀(τ, v, q, χ) ∈ V 1 × V 2 × V 3 × H 3 , ω n , τ -u n , ∇ × τ = 0, (1.4a) 1 δt u n , v + ν∇ × ω n + θω n × u n-1 + (1 -θ)ω n-1 × u n , v -p n , ∇ • v = 1 δt u n-1 , v + f n , v , (1.4b) ∇ • u n + φ n , q = 0 (1.4c) p n , χ = 0. (1.4d)
Here we used an implicit Euler discretization in time and consider a sequence of solutions at each time step indexed by an integer n ≥ 1. The real number θ ∈ [0, 1] is an arbitrary parameter. In the following we take θ = 1 2 . Let us look at a specific example in order to give some explanation on Equation (1.4): If we take

V 1 = H(curl, Ω), V 2 = H(div, Ω) and V 3 = L 2 (Ω) then we must have H 3 = {0}. Equation (1.4a) is equivalent to ω n = ∇×u n and u n ∈ H 0 (curl, Ω), (1.4b) is equivalent to u n -u n-1 δt + ν∇ × (∇ × u n ) + 1 2 (ω n × u n-1 + ω n-1 × u n ) + ∇p n = f n and p n ∈ H 1 0 , (1.4c) is equivalent to ∇ • u n = 0,
and (1.4d) is trivial here (since H 3 = {0}). This formulation is similar to the one studied by Anaya et al. [START_REF] Anaya | Analysis and Approximation of a Vorticity-Velocity-Pressure Formulation for the Oseen Equations[END_REF]. The main difference is that our formulation is studied in the framework of finite element exterior calculus and for arbitrary low order perturbations (see Equation (1.5) below). In particular the abstraction made on discrete spaces allows using any discrete subcomplex (as defined in Section 2.4). Two families are given as example in Section 2.1 but more exist on different kind of meshes (e.g. the cubical elements [START_REF] Arnold | Finite Element Exterior Calculus[END_REF]Chapter 7.7]). The construction of such families is still an active topic and the independence over the choice of discrete subcomplex is a great feature of finite element exterior calculus allowing to choose any family without modification to the proofs. Remark 1.2. Equations (1.4c) and (1.4d) might seem odd but in fact are quite natural. From the exterior calculus point of view φ and χ are harmonic forms of V 3 , in practice they can be viewed as Lagrange multipliers. For the given example

H 3 is trivial but if instead we consider V 1 = H 0 (curl, Ω), V 2 = H 0 (div, Ω) and V 3 = L 2 (Ω)
then we must have H 3 ≈ R (the 1-dimensional vector space of constant functions). Equation (1.4c) ensure that the system is onto, although here the right-hand side is null so φ n will always be zero. Equation (1.4d) ensure that p n is orthogonal to the space of harmonics 3-forms (i.e. here that Ω p n = 0).

Abstracting the linearization and time discretization scheme we simply consider two linear maps: l 3 and l 5 defined on L 2 (Ω) → L 2 (Ω) (the name of which follows the convention of Arnold [START_REF] Arnold | Finite element exterior calculus with lower-order terms[END_REF]). And we define the problem:

Given f 2 , f 3 ∈ L 2 (Ω), find (ω, u, p, φ) ∈ V 1 × V 2 × V 3 × H 3 such that ∀(τ, v, q, χ) ∈ V 1 × V 2 × V 3 × H 3 : (1.5) ω, τ -u, ∇ × τ = 0, ν∇ × ω + l 3 ω + l 5 u, v -p, ∇ • v = f 2 , v , ∇ • u, q + φ, q = f 3 , q , p, χ = 0.
We easily see that a suitable choice of l 3 and l 5 (namely (1.4). We redefine the problem in the framework of exterior calculus in (4.8). We also define a discrete counterpart to this problem in (4.14).

l 3 = (v → 1 2 v × u n-1 ) and l 5 = (v → 1 δt v + 1 2 ω n-1 × v) allows recovering
Under mild assumptions on l 3 , l 5 and Ω detailed in (4.5) we prove the wellposedness of the problem (4.8) (or equivalently (1.5)) and of its discrete counterpart (4.14). If we write (ω, u, p, φ) (resp. (ω h , u h , p h , φ h )) the solution of (4.8) (resp. (4.14)) then we derive an optimal a priori error estimate on the energy norm proportional to the approximation properties of the discrete spaces used. The result is stated in Corollary 4.17. When f 3 = 0 (which is the case for (1.4)) we show that the velocity is exactly divergence free even at the discrete level, i.e. that ∇ • u = ∇ • u h = 0 holds pointwise. We also show that the scheme is pressure-robust (see Section 5.3). This allows to derive an error estimate for the vorticity and velocity independent of the pressure in Theorem 5.5.

The remaining of this paper is divided as follows: We define the notations used in the paper and discuss some applications of the scheme in Section 2. We introduce an intermediary problem akin to a Stokes problem in Section 3. We show well-posedness and derive improved error estimates for this intermediary problem that will be useful in Section 4. Section 4 is dedicated to the analysis of problem (1.5). This is the most technical section. We derive some additional results in Section 5, and finally present a variety of numerical simulations done with our scheme in Section 6 that validate our results and give some perspectives. The exterior calculus formalism is introduced in Section 2 and heavily used in Section 3 and 4. We assume that the reader has some familiarity with exterior calculus in those two sections (3 and 4). However, no prior knowledge of exterior calculus are expected in Section 5 and 6.

2.

Setting. The exterior calculus framework allows getting a uniform vision on many objects. As such it may appear abstract and confusing at first. Therefore, we will start by giving some explicit examples of spaces to fix the ideas. We will also discuss how to deal with boundary conditions and how to use the scheme in a 2-dimensional space. Only then we will introduce the full notations and specifications of exterior calculus which will be used in the remaining of the paper.

2.1. Function spaces. Our scheme does not rely on a particular choice of discrete spaces, instead we make some assumptions (given in Section 2.4) on them and any space fulfilling these assumptions can be used. Adequate spaces are readily available on simplicial and cubic meshes (they are given e.g. in the periodic table of finite elements [START_REF] Arnold | Periodic table of the Finite Elements[END_REF]). We illustrate here an example of sensible choice. Since the discrete spaces depend on the continuous one (through boundary conditions) we begin by setting the continuous spaces for this example. Let V 1 = H 0 (curl, Ω), V 2 = H 0 (div, Ω), V 3 = L 2 (Ω) and H 3 = R ⊂ L 2 (Ω). The discrete spaces depend on a parameter r ∈ N which is a polynomial degree. Let T h be a simplicial triangulation of Ω, we choose the following discretization:

• To define the curl space V 1 h we use Nedelec's edge elements of the first kind of degree r (or

P - r Λ 1 in the periodic table), V 1 h = {ω ∈ H 0 (curl, Ω); ω |T ∈ P - r Λ 1 (T ), ∀T ∈ T h }. • To define the velocity space V 2
h we use Nedelec's face elements of the first kind of degree r (or

P - r Λ 2 ), V 2 h = {u ∈ H 0 (div, Ω); u |T ∈ P - r Λ 2 (T ), ∀T ∈ T h }. • To define the pressure space V 3
h we use discontinuous Galerkin elements of degree r -1 (or

P - r Λ 3 ), V 3 h = {p ∈ L 2 (Ω); p |T ∈ P - r Λ 3 (T ), ∀T ∈ T h }. • The space of discrete harmonic forms H 3 h = R ⊂ L 2 (Ω) is the L 2 -orthogonal complement of div(V 2 h ) in V 3 h
. We can substitute the first kind by the second, at the expense of increasing the polynomial degrees. Nedelec elements of first and second kind are respectively the 3-dimensional equivalent of Raviart-Thomas and of Brezzi-Douglas-Marini elements. The space H 3 h is just the natural way of fixing the pressure (usually only defined up to an arbitrary constant when it does not appear in boundary conditions). Figure 2.1 shows the degrees of freedom of this choice of finite elements for r = 1 and r = 2. The first element shown on the left corresponds to the space V 0 h that plays no role and will not appear in this paper.

Remark 2.1. In full generality H 3 h is simply the space of the discrete harmonic 3-forms. In this case the space of discrete harmonic 3-forms is the same as the space of continuous harmonic 3-forms. If we were to take V 1 = H 0 (curl, Ω), V 2 = H 0 (div, Ω) and V 3 = L 2 (Ω) instead then we would have H 3 = H 3 h = {0}. For this choice we can make explicit the approximation properties in terms of the size of the mesh h and of the polynomial degree r. For (ω, u, p, φ) (resp. (ω h , u h , p h , φ h )) the solution of (4.8) (resp. (4.14)) the estimate of Corollary 4.17 reads:

ω -ω h H(curl) + u -u h H(div) + p -p h L 2 + φ -φ h L 2 = O(h r ).
2.2. Boundary conditions. We give a quick review of the boundary conditions readily available. Let n be the unit outward normal to the boundary and g and h arbitrary functions in a suitable space. Equalities below are always understood on ∂Ω and any combination of the following conditions can be used:

ω × n = g, u • n = h. (2.1) ω × n = g, p = h. (2.2) u = g. (2.3) u × n = g, p = h. (2.4)
The first two conditions are uncommon since they prescribe the tangential trace of the vorticity. The third is the most commonly used since it allows enforcing the no slip condition with u • n = 0, u × n = 0.

In order to implement them, we use essential Dirichlet boundary conditions for:

• ω × n = g and take test functions in H 0 (curl).

• u • n = g and take test functions in H 0 (div). And natural conditions for the other two:

• u × n = g by adding -∂Ω (g × τ ) • nds to the left-hand side of equations (1.4).

• p = h by adding ∂Ω hv • nds to the left-hand side of equations (1.4).

Remark 2.2. These are conditions that can be easily implemented. They have the particularity of treating independently the tangential and normal components of the functions. Although they can all be imposed, only condition (2.1) and (2.4) gives a complex. The following proofs are based on the use of complexes, so we assume that one of these two conditions is chosen. When we do not use a complex we still have the well-posedness, however the rate of convergence can be impacted. The case of condition (2.3) which is used to impose the no slip condition is studied by Arnold et al. [START_REF] Arnold | Mixed finite element approximation of the vector Laplacian with Dirichlet boundary conditions[END_REF][START_REF] Chen | Convergence Analysis of Triangular MAC Schemes for Two Dimensional Stokes Equations[END_REF]. He shows that the rate of convergence on the vorticity (for the norm H(curl))

and on the pressure (for the norm L 2 ) is slightly degraded. However, for a polynomial degree of at least 2 the rate of convergence of the velocity (for the norm H(div)) is not impacted.

2.3. Two-dimensional formulation. While everything is expressed in R 3 it is perfectly possible to apply this to a 2-dimensional problem. One shall take the Raviart-Thomas's face elements (and not edge since we wish to preserve a pointwise zero divergence) or Brezzi-Douglas-Marini's face elements for

V 2 h . The meaning of ω in (1.4a) is then ω = ∇ × u = ∂ 1 u 2 -∂ 2 u 1
and the expressions appearing in (1.4b) become:

∇ × ω = ∂ 2 ω -∂ 1 ω and ω × u = ωu ⊥ = ω -u 2 u 1 .
Everything said and proved below also works in two dimensions. The independence w.r.t. dimension is a major advantage of exterior calculus formalism.

Notations.

Finally we state the full framework of exterior calculus which will be used. We make extensive use of the following notations introduced by D. Arnold [START_REF] Arnold | Finite Element Exterior Calculus[END_REF]. They are explained in greater detail in the original reference.

• d is the exterior derivative, d * its adjoint or codifferential,

• W 0 → W 1 → W 2 → W 3 is the L 2 de Rham complex of a bounded domain of R 3 .
We shall only use the last three of them,

• V 0 → V 1 → V 2 → V 3 is
a dense subcomplex on which the exterior derivative is defined (and not just densely defined). Here we have

H 1 (Ω) → H(curl, Ω) → H(div, Ω) → L 2 (Ω). • V * k is the domain of the adjoint of (d, V k ), • • is the L 2 -norm for scalar or vector valued functions, • • V is the V norm, defined by • V = • + d • ,
• Sometimes we take the norm on V ∩ V * , it is given by

• V + d * • , • V 0 h → V 1 h → V 2 h → V 3
h is a discrete subcomplex parametrized by h, • P h is the L 2 -orthogonal projection on the discrete subcomplex and in general P A is the L 2 -orthogonal projection on A. We assume that our complexes have the compactness property, which means that the inclusion V k ∩ V * k ⊂ W k is compact for each k. We also assume that there exists a cochain projection π k h :

V k → V k h , bounded for the W -norm, uniformly in h. Each space has the Hodge-decomposition W k = B k H k B * k , where B k = d(V k-1 ), B * k = d * (V * k+1 ) and H k = Z k ∩ (B k ) ⊥ , Z k being the kernel of d : V k → V k+1 .
In order to measure the approximation properties of the discrete subspaces we introduce the notation E = E k , defined by:

(2.5) ∀k, ∀σ ∈ V k , E k (σ) := inf τ ∈V k h σ -τ .
The choice of spaces introduced in Section 2.1 does indeed verify all these properties with E = O(h r ) on a dense subset (see [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF]). Each space has a Poincaré inequality (see [START_REF] Arnold | Finite Element Exterior Calculus[END_REF]Theorem 4.6]). That is to say, for each k,

∃c k p > 0 such that ∀z ∈ Z k⊥ V k , z V k ≤ c k p dz L 2 .
Moreover this also holds for the discrete subcomplex with constants bounded by c k p π k h W k (see [START_REF] Arnold | Finite Element Exterior Calculus[END_REF]Theorem 5.3]). We set c p = max k,h c k p π k h W k , so that the Poincaré inequality holds for c p regardless of the space.

When it is obvious from the context we shall drop the exponent. Hence, when applied to a 1-form v 1 ∈ V 1 we mean dv 1 := curl(v 1 ), and when applied to a 2-form v 2 ∈ V 2 we mean dv 2 := div(v 2 ). When we apply an operator (such as π h , d or d * ) to a product space, we mean to apply it to each component (i.e. d(v 1 , v 2 , v 3 ) := (dv 1 , dv 2 , 0)). When we add a suffix h such as d h instead of d we refer to the discrete counterpart of the object. We will often add a numerical suffix such as z 2 for z ∈ V 2 × V 3 , this means we take the V 2 component of z and will be clear from the context.

Usually when dealing with a primal formulation we will use the variable (u, p), and they can indeed be seen as the velocity and pressure. However, when we deal with mixed formulations we will frequently write (u 1 , u 2 , u 3 ). Here 1,2 and 3 refer to 1-form, 2-form and 3-form and have nothing to do with components in a frame of the velocity field. Specifically we will have the identification u 1 = ω, u 2 = u and u 3 = p.

The symbol A B means that there exists a constant C ∈ R * + independent of A and B (depending solely on few specified parameters) such that A ≤ CB.

3. Linear steady problem. We first study a simpler problem, analogous to a Stokes problem, which is also closely related to the Hodge-Laplace problem (see [START_REF] Arnold | Finite Element Exterior Calculus[END_REF]).

Given f = (f 2 , f 3 ) ∈ W 2 × W 3 , the problem is: Find (u, p) such that (3.1) ν∇ × (∇ × u) + ∇p = f 2 on Ω, ∇ • u = f 3 on Ω.
3.1. Primal formulation. We give the primal formulation in order to explain more clearly what we seek in the mixed formulation, and because we shall need this operator in Section 4. Since we have no way to reach the harmonic part of f 2 , we must include a second harmonic space (this time of 2-forms).

Definition 3.1. Let D 0 := {(u, p) ∈ (V 2 ∩ V * 2 ) × V * 3 | d * u ∈ V 1 } , f ∈ W 2 × W 3 and (3.2) L 0 := νdd * -d * d 0 = ν∇ × ∇× ∇ ∇• 0 .
The problem is to find (u, p) ∈ P H ⊥ D 0 , such that

(3.3) ∀(v, q) ∈ P H ⊥ D 0 , L 0 (u, p), (v, q) = f, (v, q) .
Remark 3.2. Since d : V 3 → 0, the Hodge decomposition on W 3 reads W 3 = B 3 H 3 . We also have the equality between the V -norm and the W -norm.

Continuous

Well-posedness. Now we introduce the mixed formulation and switch to the exterior calculus notation. Recall that u 2 (resp. u 3 ) defined below corresponds to u (resp. p) in (3.1). The mixed formulation is characterized by the introduction of the auxiliary variable u 1 corresponding to ∇ × u 2 or d * u 2 . The problem reads:

Given (f 2 , f 3 ) ∈ W 2 × W 3 , find (u 1 , u 2 , u 3 , φ 2 , φ 3 ) ∈ V 1 × V 2 × V 3 × H 2 × H 3 such that ∀(v 1 , v 2 , v 3 , χ 2 , χ 3 ) ∈ V 1 × V 2 × V 3 × H 2 × H 3 , (3.4) u 1 , v 1 -u 2 , dv 1 = 0, ν du 1 , v 2 -u 3 , dv 2 + φ 2 , v 2 = f 2 , v 2 , du 2 , v 3 + φ 3 , v 3 = f 3 , v 3 , u 2 , χ 2 + u 3 , χ 3 = 0.
The associated bilinear form is noted B 0 . Remark 3.3. There are two important reasons to remove d * from the formulation. They become clear when we will consider the discrete counterpart of the formulation. The first one is that d * and d * h are barely related, and it is much harder to analyze the difference between them than between d and d h . The second one is because while d h is a local operator, d * h is in general a global operator. This greatly deteriorates the sparsity pattern of the system.

Lemma 3.4. There is α > 0, such that ∀(u 1 , u 2 , u 3 , φ 2 , φ 3 ) ∈ V 1 ×V 2 ×V 3 ×H 2 ×H 3 we have: sup (v1,v2,v3,χ2,χ3) B 0 ((u 1 , u 2 , u 3 , φ 2 , φ 3 ), (v 1 , v 2 , v 3 , χ 2 , χ 3 )) (v 1 , v 2 , v 3 , χ 2 , χ 3 ) V ≥ α (u 1 , u 2 , u 3 , φ 2 , φ 3 ) V .
Proof. For any (u

1 , u 2 , u 3 , φ 2 , φ 3 ) ∈ V 1 × V 2 × V 3 × H 2 × H 3 , let ρ 1 ∈ B * 1 ∩ V 1 be such that dρ 1 = P B u 2 , ρ 2 ∈ B * 2 ∩ V 2 be such that dρ 2 = P B u 3 + du 2 .
And take

v 1 = νu 1 -ν c 2 p ρ 1 , v 2 = P B u 2 + du 1 -ρ 2 + φ 2 , v 3 = du 2 -P B u 3 + φ 3 , χ 2 = P H u 2 , χ 3 = P H u 3 .
The Poincaré inequality gives:

ρ 1 V 1 ≤c p P B u 2 , (3.5) ρ 2 V 2 ≤c p P B u 3 + c p du 2 . (3.6)
And we easily see that:

v 1 V 1 + v 2 V 2 + v 3 V 3 + χ 2 + χ 3 u 1 V 1 + u 2 V 2 + u 3 V 3 + φ 2 + φ 3 ,
where the hidden constant only depends on ν and c p . Using orthogonality of the Hodge decomposition we get:

(3.7)

B 0 ((u 1 , u 2 , u 3 , φ 2 , φ 3 ), (v 1 , v 2 , v 3 , χ 2 , χ 3 )) = u 1 , νu 1 - ν c 2 p ρ 1 -u 2 , d(νu 1 - ν c 2 p ρ 1 ) + ν du 1 , P B u 2 + du 1 -ρ 2 + φ 2 -u 3 , d(P B u 2 + du 1 -ρ 2 + φ 2 ) + du 2 , du 2 -P B u 3 + φ 3 + u 2 , P H u 2 + u 3 , P H u 3 + φ 2 , P B u 2 + du 1 -ρ 2 + φ 2 + φ 3 , du 2 -P B u 3 + φ 3 = ν u 1 , u 1 - ν c 2 p u 1 , ρ 1 -ν u 2 , du 1 + ν c 2 p u 2 , P B u 2 + ν du 1 , u 2 + ν du 1 , du 1 + u 3 , P B u 3 + u 3 , du 2 + du 2 , du 2 -du 2 , u 3 + u 2 , P H u 2 + u 3 , P H u 3 + φ 2 , φ 2 + φ 3 , φ 3 = ν u 1 2 - ν c 2 p u 1 , ρ 1 + ν c 2 p P B u 2 2 + ν du 1 2 + P B u 3 2 + du 2 2 + P H u 2 2 + P H u 3 2 + φ 2 2 + φ 3 2 . Now since u 1 , ρ 1 ≤ u 1 ρ 1 ≤ c 2 p 2 u 1 2 + 1 2c 2 p ρ 1 2
, using (3.5) and injecting in the expression (3.7) we get:

B 0 ((u 1 ,u 2 , u 3 , φ 2 , φ 3 ), (v 1 , v 2 , v 3 , χ 2 , χ 3 )) ≥ ν 2 u 1 2 + ν 2c 2 p P B u 2 2 + ν du 1 2 + u 3 2 + du 2 2 + P H u 2 2 + P H u 3 2 + φ 2 2 + φ 3 2 .
Finally, since du 3 = 0, dφ 2 = dφ 3 = 0 and P B * u 2 ≤ c p dP B * u 2 = du 2 we have:

B 0 ((u 1 , u 2 , u 3 , φ 2 , φ 3 ), (v 1 , v 2 , v 3 , χ 2 , χ 3 )) (u 1 , u 2 , u 3 , φ 2 , φ 3 ) 2 V ,
where the hidden constant depends only on c p and ν.

Lemma 3.5. Given any

(v 1 , v 2 , v 3 , χ 2 , χ 3 ) ∈ V 1 × V 2 × V 3 × H 2 × H 3 there is (u 1 , u 2 , u 3 , φ 2 , φ 3 ) ∈ V 1 × V 2 × V 3 × H 2 × H 3 such that B 0 ((u 1 , u 2 , u 3 , φ 2 , φ 3 ), (v 1 , v 2 , v 3 , χ 2 , χ 3 )) > 0 . Proof. If P B v 2 = 0 take u 1 ∈ B * 1 such that du 1 = P B v 2 , φ 2 = 0, u 3 = φ 3 = 0. Then if u 1 , v 1 = 0 take u 2 = 0 else take u 2 = dv 1 u1,v1 dv1,dv1 (dv 1 = 0 since P B * v 1 = 0 as u 1 , v 1 = 0). If P B v 2 = 0, simply take u 1 = v 1 , u 2 such that du 2 = P B v 3 , P B u 2 = -dv 1 , P H u 2 = χ 2 (
this is possible by the Hodge decomposition),

u 3 such that P B u 3 = -dv 2 , P H u 3 = χ 3 , φ 2 = P H v 2 and φ 3 = P H v 3 .
Lemma 3.4 together with Lemma 3.5 give the conditions to apply the Babuška-Lax-Milgram theorem. This proves the continuous well-posedness. Moreover, we have:

(3.8) u 1 V 1 + u 2 V 2 + u 3 V 3 + φ 2 V 2 + φ 3 V 3 ≤ c f
where c depends only on c p and ν.

Discrete well-posedness.

The discrete problem reads:

Given (f 2 , f 3 ) ∈ W 2 × W 3 , find (u 1h , u 2h , u 3h , φ 2h , φ 3h ) ∈ V 1 h × V 2 h × V 3 h × H 2 h × H 3 h such that ∀(v 1h , v 2h , v 3h , χ 2h , χ 3h ) ∈ V 1 h × V 2 h × V 3 h × H 2 h × H 3 h , (3.9) B 0 ((u 1h , u 2h , u 3h , φ 2h , φ 3h ), (v 1h , v 2h , v 3h , χ 2h , χ 3h )) = f 2 , v 2h + f 3 , v 3h .
Since we have a discrete Poincaré inequality (see [START_REF] Arnold | Finite Element Exterior Calculus[END_REF]) and use a subcomplex we can apply exactly the same proof as in the continuous case, substituting continuous spaces for their discrete counterparts.

Remark 3.6. The only differential operator used in the mixed formulation is d and its discrete counterpart is merely its restriction to the discrete space. However since the discrete space of harmonic forms H h is not necessarily a subspace of H. When H 2 or H 3 is not zero, the method may be non-conforming.

3.4. Error estimate. We will only derive a basic error estimate for a global norm here. First we define:

(3.10) µ := max k∈{2,3} sup r∈H k , r =1 (I -π k h )r . Theorem 3.7. Given (f 2 , f 3 ) ∈ W 2 × W 3 , let (u 1 , u 2 , u 3 , φ 2 , φ 3
) be the solution of the continuous problem (3.4) and (u 1h , u 2h , u 3h , φ 2h , φ 3h ) the solution of the discrete problem (3.9), then for E given by (2.5) it holds:

(u 1 -u 1h , u 2 -u 2h , u 3 -u 3h ) V + (φ 2 -φ 2h , φ 3 -φ 3h ) inf v1 ∈V 1 h u 1 -v 1 V 1 + inf v2 ∈V 2 h u 2 -v 2 V 2 + E(u 3 ) + E(φ 2 ) + E(φ 3 ) + µ(E(P B u 2 ) + E(P B u 3 )). Proof. For all (v 1h , v 2h , v 3h , χ 2h , χ 3h ) ∈ V 1 h × V 2 h × V 3 h × H 2 h × H 3 h we have: B 0 ((u 1 , u 2 , u 3 , φ 2 , φ 3 ), (v 1h ,v 2h , v 3h , χ 2h , χ 3h )) = f 2 , v 2h + f 3 , v 3h + u 2 , χ 2h + u 3 , χ 3h .
Let (v 1 , v 2 , v 3 , χ 2 , χ 3 ) be the V -orthogonal projection of (u 1 , u 2 , u 3 , φ 2 , φ 3 ) into their respective discrete spaces. By the continuity of B 0 it holds,

∀(v 1h , v 2h , v 3h , χ 2h , χ 3h ) ∈ V 1 h × V 2 h × V 3 h × H 2 h × H 3 h B 0 ((u 1h -v 1 , u 2h -v 2 , u 3h -v 3 , φ 2h -χ 2 , φ 3h -χ 3 ), (v 1h , v 2h , v 3h , χ 2h , χ 3h )) = B 0 ((u 1 -v 1 , u 2 -v 2 , u 3 -v 3 , φ 2 -χ 2 , φ 3 -χ 3 ), (v 1h , v 2h , v 3h , χ 2h , χ 3h )) -u 2 , χ 2h -u 3 , χ 3h = B 0 ((u 1 -v 1 , u 2 -v 2 , u 3 -v 3 , φ 2 -χ 2 , φ 3 -χ 3 ), (v 1h , v 2h , v 3h , χ 2h , χ 3h )) -P H h u 2 , χ 2h -P H h u 3 , χ 3h ( u 1 -v 1 V 1 + u 2 -v 2 V 2 + u 3 -v 3 + φ 2 -χ 2 + φ 3 -χ 3 + P H h u 2 V 2 + P H h u 3 V 3 )( v 1h V 1 + v 2h V 2 + v 3h + χ 2h + χ 3h ),
where the hidden constant only depends on ν.

From the discrete inf-sup condition we have

( u 1h -v 1 V 1 + u 2h -v 2 V 2 + u 3h -v 3 + φ 2h -χ 2 + φ 3h -χ 3 ) ( u 1 -v 1 V 1 + u 2 -v 2 V 2 + u 3 -v 3 + φ 2 -χ 2 + φ 3 -χ 3 + P H h u 2 V 2 + P H h u 3 ),
where the hidden constant only depends on ν and the discrete constant of Poincaré.

The theorem follows from

P H h u i V µE(P B u i ), ∀i ∈ {2, 3}.
This is proved in Theorem 5.2 of [START_REF] Arnold | Finite Element Exterior Calculus[END_REF].

Remark 3.8. E(φ i ) is understood as viewing φ i as an element of

V i ⊃ H i , E(φ i ) = inf q∈V i h φ i -q .
4. Linearized problem. Starting from the linear steady problem (3.4), we add some terms of lower order since this will allow us to construct our scheme. Again we recall the correspondence with the names used for variables in the introduction: (u 1 , u 2 , u 3 , u p ) ≡ (ω, u, p, φ) in (1.5). To keep the notations bearable we shall also write u B = P B u, u H = P H u and so on.

We consider two linear maps l 3 : W 1 → W 2 and l 5 : W 2 → W 2 (we chose these names to match those used by D. Arnold [START_REF] Arnold | Finite element exterior calculus with lower-order terms[END_REF]). Define

D = {(u, p) ∈ (V 2 ∩ V * 2 ) × (V * 3 ∩ H 3⊥ )| d * u ∈ V 1 } , W = W 2 × (W 3 ∩ H 3⊥ ), (4.1) L := (νd + l 3 )d * + l 5 -d * d 0 .
We consider the primal problem: Given f ∈ W , to find (u, p) ∈ D such that (4.2) L(u, p) = f.

We also define the dual operator

L on D = {(u, p) ∈ (V 2 ∩ V * 2 ) × (V * 3 ∩ H 3⊥ )| (νd * + l * 3 )u ∈ V 1 } by: (4.3) L := d(νd * + l * 3 ) + l * 5 d * -d 0
As an intermediary step, we wish to extend L on a larger domain and introduce the following notations: (4.4)

L λ : (V 2 ∩ V * 2 ) × (V * 3 ∩ H 3⊥ ) → (V * 2 × (W 3 ∩ H 3⊥ )) , L λ (u, p)(v, q) := νd * u, d * v + l 3 d * u + l 5 u -λd * p, v + λdu, q , L λ : (V 2 ∩ V * 2 ) × (V * 3 ∩ H 3⊥ ) → (V * 2 × (W 3 ∩ H 3⊥ )) , L λ (u, p)(v, q) := νd * u + l * 3 u, d * v + l * 5 u + λd * p, v + -λdu, q .
Here λ is a positive parameter introduced to simplify the proof of Theorem 4.1. In the Theorem 4.1 we shall see that they are almost always isomorphisms. We define the solution operator K = (L 1 ) -1 , and we assume that

(4.5) d * (L 1 -1 ) 2 (W ) ⊂ V 1 , (νd * + l * 3 )(K) 2 (W ) ⊂ V 1
where (L 1 -1 ) 2 and (K) 2 are the projections on the first component of the product space taken after the operators. Moreover, we assume that dd * (K) 2 W →W 2 and dl * 3 (K) 2 W →W 2 are bounded. We show in Section 5.1 that these assumptions are satisfied when l 3 and l 5 are those used in our scheme.

The proof follows the same outline as [START_REF] Arnold | Finite element exterior calculus with lower-order terms[END_REF]. First we prove that the continuous primal formulation gives an isomorphism, then we prove that the continuous mixed formulation is well-posed. Lastly we prove the well-posedness of the discrete mixed formulation and give an estimation of the error in energy norm.

Continuous primal formulation.

Theorem 4.1. L λ + µ •, • is a bounded isomorphism for all µ ∈ C except for a countable subset.

Proof. Let c = max( l 3 , l 5 ), for (u, p) ∈ (V 2 ∩ V * 2 ) × (V * 3 ∩ H 3⊥ ) take v B = u B , v H = u H , v B * = -d * p, q = du. We have: (4.6) (µ •, • + L λ )(u, p)(v, q) = ν d * u, d * u + λ d * p, d * p + λ du, du + l 3 d * u, u B -d * p + u H + l 5 u, u B -d * p + u H + µ u, u B + u H -µ u, d * p + µ p, du = ν d * u 2 + λ d * p|| 2 + λ du 2 + µ u B + u H 2 + l 3 d * u, u B -d * p + u H + l 5 u, u B -d * p + u H .
With the Cauchy-Schwarz inequality and the identity ab ≤ 2 -1 (a 2 + b 2 ) we bound the last line from (4.6):

| l 3 d * u,u B -d * p + u H | + | l 5 u, u B -d * p + u H | ≤ c 2 2ν ( u B + u H 2 + d * p 2 ) + ν 2 d * u 2 + c 2 2 ( u B + u H 2 + d * p 2 ) + 1 2 u 2 .
Since u B * 2 ≤ c 2 p du 2 and u 2 = u B * 2 + u B + u H 2 we have:

(µ •, • + L λ )(u, p)(v, q) ≥ (µ - c 2 2ν - c 2 2 - 1 2 ) u B + u H 2 + ν 2 d * u 2 + (λ - c 2 2ν - c 2 2 ) d * p|| 2 + (λ - c 2 p 2 ) du 2 .
Using once again the Poincaré inequality to bound u B * by du and p B by d * p (on the dual complex), and since p B = p (as p ∈ V 3 ∩ H 3⊥ ) we have for λ and µ large enough (solely depending on c, ν and c p ):

(µ •, • + L λ )(u, p)(v, q) u 2 + du 2 + d * u 2 + p 2 + d * p 2 (u, p) 2 V 2 ∩V * 2 ×V * 3 . Clearly (v, q) V * 2 ×W 3 ≤ (u, p) V 2 ∩V * 2 ×V * 3 as d * v = d * u B ≤ d * u and (µ •, • + L λ ) is continuous as a bilinear form from (V 2 ∩ V * 2 × V * 3 ) × (V * 2 × (W 3 ∩ H 3⊥ ))
. The only thing left to show in order to use the Babuška-Lax-Milgram theorem is the second condition. For any (v, q) = 0 we must find (u, p) such that (µ •, • + L λ )(u, p)(v, q) > 0. We can take λ and µ such that λ 2 = µν. We consider two cases: [START_REF] Arnold | Finite element exterior calculus with lower-order terms[END_REF]). Take p = -ν/λdu, then since λ = µν/λ it holds:

When v = 0, we can find u ∈ {w ∈ V 2 ∩ V * 2 |d * w ∈ V 1 , dw ∈ V * 3 } such that ((νd + l 3 )d * + l 5 + µ + νd * d)u = v (see
(µ •, • + L λ )(u, p)(v, q) = νdd * u, v + (l 3 d * + l 5 + µ)u, v -λ d * p, v + λ du, q + µ p, q = ((νd + l 3 )d * + l 5 + µ + νd * d)u, v + λ du, q - µν λ du, q = v, v .
When v = 0 take u solution of the Hodge-Dirac problem (see [START_REF] Leopardi | The Abstract Hodge-Dirac Operator and Its Stable Discretization[END_REF]): du = q, d * u = 0. Then (µ •, • + L λ )(u, 0)(0, q) = λ du, q = λ q, q > 0.

This shows that (µ

•, • +L λ ) is a bounded isomorphism from (V 2 ∩V * 2 )×(V * 3 ∩H 3⊥ ) to (V * 2 ×(W 3 ∩H 3⊥ )) . Since I : (V 2 ∩V * 2 )×(V * 3 ∩H ⊥ 3 ) → (V * 2 ×(W 3 ∩H 3⊥
)) is compact by the compactness property, I(µ •, • + L λ ) -1 is also compact. Since the spectrum of a compact operator is at most countable, we have that Id + ηI(µ •, • + L λ ) -1 has a bounded inverse for all η ∈ C except for a countable subset. Therefore, by composing to the right with µ •, • + L λ we get that L λ + (µ + η)I has almost always a bounded inverse.

Hence, up to an arbitrary small perturbation, L λ is a bounded isomorphism from

(V 2 ∩ V * 2 ) × (V * 3 ∩ H ⊥ 3 ) to (V * 2 × (W 3 ∩ H 3⊥ )) . Remark 4.2.
We could have left H 3 in the domain and the proof above would still work. However, in this case, L λ would never have been an isomorphism since its image cannot reach H 3 .

We have the same result for the dual problem. Lemma 4.3. For all µ ∈ C except for a countable subset,

L λ + µ •, • is a bounded isomorphism. Lemma 4.8. For any v ∈ V , there is u ∈ V such that B(u, v) > 0. Proof. Given v = 0 ∈ V , if v 2 = 0, v 3 =
0 and v p = 0 take u = (v 1 , 0, 0, 0) and

B(u, v) = v 1 , v 1 > 0. Else take (u 2 , u 3 ) = L 1 -1 (v 2 , P H ⊥ v 3 ) + (0, v p ), u p = P H v 3 and u 1 = d * u 2 (u 1 ∈ V 1 by assumption (4.5)) then B(u, v) = v 2 , v 2 + v 3 , v 3 + q, q > 0.
4.3. Discrete well-posedness. We introduce the notation

V h = V 1 h × V 2 h × V 3 h × H 3 h .
The discrete variational problem is the same as the continuous, substituting V by V h . Hence we shall still use the notation B, this time as a function from

V h × V h to R. So that the discrete problem is: Given (f 2 , f 3 ) ∈ W , find u h = (u 1h , u 2h , u 3h , u ph ) ∈ V 1 h × V 2 h × V 3 h × H 3 h such that ∀v h = (v 1h , v 2h , v 3h , v ph ) ∈ V 1 h × V 2 h × V 3 h × H 3 h , (4.14) B(u h , v h ) = f 2 , v 2h + f 3 , v 3h .
Considering the dual problem to the unperturbed problem with ν = 1. We have:

D 0 = {(u, p) ∈ (V 2 ∩ V * 2 ) × V 3 * | d * u ∈ V 1 } and L 0 (u, p) = (dd * u + d * p, -du).
Let K 0 be the solution operator of the dual problem L 0 . We have K 0 = (L 0 ) -1 when L 0 is viewed as an isomorphism from P H ⊥ D 0 to P H ⊥ (W 2 × W 3 ) and K 0 is extended by 0 on H. Explicitly we have the decomposition:

∀(f 2 , f 3 ) ∈ W 2 × W 3 , (4.15) (f 2 , f 3 ) = (dd * (K 0 ) 2 (f 2 , f 3 ) + d * (K 0 ) 3 (f 2 , f 3 ), -d(K 0 ) 2 (f 2 , f 3 )) + (P H f 2 , P H f 3 )
and a similar expression for their discrete counterparts. Therefore, ∀(z 2 , z 3 ) ∈ D 0 , (P

H ⊥ z 2 , P H ⊥ z 3 ) = L 0 K 0 (z 2 , z 3 ) = K 0 L 0 (z 2 , z 3 ).
As mentioned before this problem is closely related to the one studied by D. Arnold [START_REF] Arnold | Finite element exterior calculus with lower-order terms[END_REF]. Since the mixed variable part is almost unchanged we shall use the generalized canonical projection Π h (see [START_REF] Arnold | Finite element exterior calculus with lower-order terms[END_REF]) and we state its properties below. Lemma 4.9. Under the condition of [7, Theorem 5.1]:

• Π h is a projection uniformly bounded in the V-norm.

• dΠ h = P B h d. • ∀w ∈ V k , (I -Π h )w (I -π h )w + η 0 dw . • ∀w, v ∈ V k , | (I -Π h )w, v | ( (I -π h )w + η 0 dw )( (I -π h )v + η 0 dv ) + α 0 dw dv .
where η 0 , α 0 → 0 when h → 0, they are given, along the proof in the reference [START_REF] Arnold | Finite element exterior calculus with lower-order terms[END_REF].

Definition 4.10. We shall use the following notations in this section:

δ 0 = (I -π h )K 0 , µ 0 = (I -π h )P H , η 0 = max{ (I -π h )dK 0 , (I -π h )d * (K 0 ) 2 }, α 0 = η 0 (1 + η 0 ) + µ 0 + δ 0 + µ 0 δ 0 + η 0 , η = max{δ 0 , µ 0 , η 0 , (I -π h )l * 3 (K) 2 , (I -π h )dl * 3 (K) 2 }. Lemma 4.11.
We have:

K 0 -K 0h P h α 0 , dK 0 -dK 0h P h + d * (K 0 ) 2 -d * h (K 0h ) 2 P h η 0 .
Proof. The idea is for (f 2 , f 3 ) ∈ (W 2 × W 3 ) to apply the error estimate of Theorem 3.7 for (u 2 , u

3 ) = K 0 (f 2 , f 3 ), (φ 2 , φ 3 ) = P H (f 2 , f 3 ), u 1 = d * u 2 , (u 2h , u 3h ) = K 0h P h (f 2 , f 3 ), (φ 2h , φ 3h ) = P H h P h (f 2 , f 3 ), u 1h = d * h u 2h .
Unfortunately we cannot conclude with the crude estimate of Theorem 3.7 because of the error on u 1 . We need improved estimates that give

u 2 -u 2h + u 3 -u 3h (1 + µ 0 )E(u 2 ) + E(u 3 ) + η 0 E(u 1 ) (η 2 0 + δ 0 + η 0 )E(du 1 ) + η 0 E(φ 2 ) + E(du 2 ), du 2 -du 2h + u 1 -u 1h E(du 2 ) + E(u 1 ) + η 0 E(du 1 ). We conclude since E(u 2 ) + E(u 3 ) ≤ δ 0 (f 2 , f 3 ) , E(u 1 ) + E(du 2 ) ≤ η 0 (f 2 , f 3 ) and E(du 1 ) + E(φ 2 ) ≤ (f 2 , f 3
) , the last coming from du 1 = P B f 2 . These proofs are lengthy, technical and mostly follow those in theorem 3.11 of [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF].

Lemma 4.12. For f ∈ H ⊥ we have

P H h f µ 0 f
Proof. We recall the mixed formulation for the Hodge-Laplacian problem. The bilinear form is given by: B((σ, u, p), (τ, v, q)) = σ, τ -u, dτ + dσ, v + du, dv + p, v + u, q .

In the continuous case the bilinear form acts on (

V 1 ×V 2 ×H 2 ) 2 and on (V 1 h ×V 2 h ×H 2 h ) 2 in the discrete case. Let (σ, u, p) ∈ V 1 × V 2 × H 2 be such that ∀(τ, v, q) ∈ V 1 × V 2 × H 2 , B((σ, u, p), (τ, v, q)) = (f, v), and let (σ h , u h , p h ) ∈ V 1 h × V 2 h × H 2 h be such that ∀(τ, v, q) ∈ V 1 h × V 2 h × H 2 h , B((σ h , u h , p h ), (τ, v, q)) = (f, v).
Then [4, Theorem 5.6 p. 62] gives the error estimate p -p h E(p) + µ 0 E(dσ). We have P H f = p = 0 and

P H h f = p h thus P H h f = p -p h 0 + µ 0 E(dσ) µ 0 f
since E(p) = 0 and by the well-posedness of the Hodge-Laplacian problem for the last inequality.

Theorem 4.13. For z = (z 2 , z 3 ) ∈ D 0 , let z h = (z 2h , z 3h ) = K 0h P h L 0 z +P H h P H z, we have:

z -z h α 0 L 0 z , d(z -z h ) + d * z 2 -d * h z 2h η 0 L 0 z , P h (dd * z 2 + d * z 3 ) -(dd * h z 2h + d * h z 3h ) ≤ µ 0 L 0 z . Proof.
The same proof as [7, Theorem 5.2] works. It is duplicated here since the proof is short and demonstrates that our change in the definition of L 0 does not interfere.

z -z h = (P H ⊥ z -P H ⊥ h z h ) + (P H z -P H h z h ) = (K 0 -K 0h P h )L 0 z + (I -P H h )P H z.
As P H h P H = P Z h P H and π h Z ⊂ Z h we have:

(I -P H h )P H z ≤ (I -π h )P H z ≤ µ 0 z α 0 L 0 z .
We then get the expected result by Lemma 4.11. The second part follows directly from Lemma 4.11. Finally for the last estimate, (4.15) gives

dd * h z 2h + d * h z 3h = (L 0h ) 2 z h = (L 0h K 0h P h ) 2 L 0 z = (P B h + P B * h )P h (L 0 ) 2 z so: P h (dd * z 2 + d * z 3 ) -(dd * h z 2h + d * h z 3h ) = (I -(P B h + P B * h ))P h (L 0 ) 2 z = P H h (L 0 ) 2 z .
And we conclude with Lemma 4.12.

Given u ∈ V h , define g = (u 2 , P H ⊥ u 3 ), z = Kg, ξ = -(d * + l * 3 )z 2 and z = (ξ, z, P H u 3 ).

Theorem 4.14. There is

z h ∈ V h such that ∀ω ∈ V h , z h V z V uniformly in h and |B(ω, z -z h )| ≤ h ω V u , where h → 0 as h → 0. Proof. Take z h = K 0h P h L 0 z + P H h P H z, ξ h = -d * h z 2h -Π h l * 3 z 2 .
By Theorem 4.13 we have:

(4.16) z -z h α 0 g , d(z -z h ) η 0 g , d * z 2 -d * h z 2h η 0 g , ξ -ξ h ≤ d * z 2 -d * h z 2h + (I -Π h )l * 3 z 2 .
Using Lemma 4.9 and the boundedness of d l * 3 K we get

(I -Π h )l * 3 z 2 (I -π h )l * 3 z 2 + η 0 dl * 3 z 2 (η + η 0 ) g . Finally, since ∀ω 2 ∈ V 2 h ⊂ V 2 dω 2 , z 3 = ω 2 , d * z 3 , dω 2 , z 3h = ω 2 , d * h z 3h ,
we have ∀ω ∈ V h :

|B(ω,(ξ -ξ h , z 2 -z 2h , z 3 -z 3h , P H u 3 -P H h P H u 3 ))| = | ω 1 , ξ -ξ h -ω 2 , d(ξ -ξ h ) + (νd + l 3 )ω 1 , z 2 -z 2h -ω 3 , d(z 2 -z 2h ) + l 5 ω 2 , z 2 -z 2h + dω 2 , z 3 -z 3h + ω p , z 3 -z 3h + ω 3 , P H u 3 -P H h P H u 3 | ω V ( ξ -ξ h + 2 z 2 -z 2h + d(z 2 -z 2h ) + z 3 -z 3h + (I -P H h )P H u 3 ) + | ω 2 , P h d * z 3 -d * h z 3h -P h d(ξ -ξ h ) | ω V [(η + 2η 0 + 2α 0 + η 0 + µ 0 ) u + P h d * z 3 -d * h z 3h -P h d(ξ -ξ h ) ] .
Since η, η 0 , η 0 , µ 0 and α 0 all converge toward 0 when h → 0, the only thing left to prove is that

P h d * z 3 -d * h z 3h -P h d(ξ -ξ h )
u where → 0 when h → 0. We have

-d(ξ -ξ h ) = d(d * + l * 3 )z 2 -dd * h z 2h -dΠ h l * 3 z 2 = dd * z 2 -dd * h z 2h + d(I -Π h )l * 3 z 2 .
Theorem 4.13 gives:

P h d * z 3 -d * h z 3h + P h dd * z 2 -dd * h z 2h µ 0 g .
And we conclude with Lemma 4.9 since we can find a bounded cochain projection π h such that π h d = P B h d (see [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF]Theorem 3.7]) so

d(I -Π h )l * 3 z 2 = (I -P B h )dl * 3 z 2 (I -π h )dl * 3 z 2 ≤ η g .
Lemma 4.15. For all u ∈ V h and z ∈ V defined in Theorem 4.14, there exists a constant c independent of h and σ ∈ V h such that σ V u V and B(u, z + σ) ≥ c u 2 V . Proof. Starting from Lemma 4.6, we construct σ in the same way as we did in Lemma 4.7 in the continuous case. We must simply add B(ω, (0, 0, -P H h z 3 , 0)) = -ω p , z 3 to correct the harmonic part.

Theorem 4.16. There are two positive constants h 0 and C 0 such that for all h ∈ (0, h 0 ], there exists a unique u ∈ V h such that ∀v ∈ V h , B(u, v) = (f , v). Moreover we have u V ≤ C 0 f . Proof. For u ∈ V h and z, z h defined in Theorem 4.14, Lemma 4.15 gives σ ∈ V h with σ V u V and a constant c independent of h such that:

|B(u, z + σ)| ≥ c u 2 V .
By Theorem 4.14 we have for a constant b independent of h.

B(u, z h -z) ≤ h b u 2 V .
Therefore:

|B(u, σ + z h )| = |B(u, z + σ) + B(u, z h -z)| ≥ |B(u, z + σ)| -|B(u, z h -z)| ≥ c u 2 V -h b u 2 V ≥ (c -h b) u 2 V .
Since h → 0 as h → 0 we can find h 0 such that ∀h ∈ (0, h 0 ], ch b > c -h0 b > 0. By Theorem 4.14 and by the expression of σ we find:

σ + z h V u V + z V u V .
This ends the proof since V h has finite dimension.

Corollary 4.17. If assumption (4.5) holds then for h ≤ h 0 given by Theorem 4.16, and for u (resp. u h ) the solution of the continuous problem (4.8) (resp. of the discrete problem (4.14)) it holds:

u -u h V inf v1∈V 1 h u 1 -v 1 V 1 + inf v2∈V 2 h u 2 -v 2 V 2 + E(u 3 ) + E(u p ) + µ 0 E(P B u 3 ).
Proof. The proof is the same as in Theorem 3.7

Remark 4.18. The hidden constant of Corollary 4.17 depends on l 5 which, in the case of problem (1.4), blows up when δt → 0. A more subtle analysis is required to explicit the dependency of the error on δt. For a single time step in the setting of Corollary 4.17 the error will actually decrease when δt → 0.

5. Conserved quantities. Lastly we prove that our scheme does indeed verify properties mentioned in the introduction as well as the regularity assumption (4.5).

Regularity assumptions.

Problem (1.4) is a special case of Problem (4.8) taking suitable l 3 and l 5 . We prove below that assumptions (4.5) are valid if u n-1 ∈ H 2 (Ω) and if the domain is smooth enough to have H 0 (curl, Ω) ∩ H(div, Ω) ⊂ H 1 (Ω), as discussed in [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF]Chapter 3.2]. We use the notation H 1 (Ω) both for scalar and for vector fields. First we need the following lemma:

Lemma 5.1. If B ∈ H 1 (Ω) and A ∈ H 2 (Ω) then A × B ∈ H 1 (Ω) for a smooth enough domain Ω of R 3 (or R 2 ).
Proof. We have H 1 (Ω) ⊂ L 4 (Ω) and H 2 (Ω) ⊂ C 0 (Ω) by Sobolev Embedding theorems thus ∀i, j, k ∈ {x, y, z}, ∂ i A j ∈ L 4 (Ω), ∂ i B j ∈ L 2 (Ω) and terms of the form A i ∂ j B k are the product of a bounded function with a function in L 2 (Ω) and those of the form B i ∂ j A k are the product of two functions in L 4 (Ω).

Going back to problem (1.4), we take

V 1 = H(curl, Ω), V 2 = H(div, Ω), V 3 = L 2 (Ω) and H 3 = 0. Regarding assumption (4.5), K is a bounded isomorphism from L 2 × L 2 to (H(div) ∩ H 0 (curl)) × H 1 0 ⊂ H 1 × H 1 .
If we assume u ∈ H 2 and v ∈ H 1 then using the scalar triple product we have ∀σ ∈ H(curl):

2 l 3 σ, v = (σ × u n-1 ) • v = (u n-1 × v) • σ = 2 σ, l * 3 v . Thus l * 3 v = 1 2 u n-1 × v on H 1 , l * 3 maps H 1 into itself by Lemma 5.1, and
l * 3 (K) 2 (W ) ⊂ H 1 ⊂ H(curl) = V 1 . Thus we have (∇×)l * 3 (K) 2 W →L 2 (K) 2 W →H 1 and by the boundedness of K L 2 ×L 2 →H 1 ×H 1 we get the boundedness of dl * 3 (K) 2 W →L 2 . Finally we have νdd * (K) 2 = I -dl * 3 (K) 2 + l * 5 (K) 2 + d * (K) 3 as distributions.
From the L 2 → L 2 boundedness of the right-hand side we get both νd * (K) 2 (W ) ⊂ H(curl) = V 1 and the boundedness of dd * (K) 2 W →L 2 .

The same argument applied to L

-1 shows that d * (L -1 ) 2 (W ) ⊂ V 1 .
Hence (4.5) is fulfilled.

Remark 5.2. Assuming u n-1 ∈ H 2 (Ω) is very mild as any solution u of (4.2) must have ∇ × u ∈ H(curl, Ω), ∇ • u = 0 thus ∆u ∈ L 2 (Ω). Hence by elliptic regularity for Ω smooth enough and if u satisfies appropriate boundary conditions then u ∈ H 2 . 5.2. Pointwise vanishing divergence. This is a simple fact that follows from the use of a discrete subcomplex. The operator d h = div used in the formulation is simply the restriction of the continuous operator to the discrete space V 2 h and its image is contained in the discrete space V 3 h . This holds even for condition (2.2) and (2.3) which does not give a complex structure. From (1.4) we have ∀q

h ∈ V 3 h , ∇ • u h + φ h , q h = 0 with ∇ • u h ⊥ φ h and ∇ • u h ∈ V 3 h by construction. Therefore, taking q h = ∇ • u h we have ∇ • u h + φ h , ∇ • u h = ∇ • u h 2 = 0.
5.3. Pressure-robustness. A scheme is called pressure-robust ( [START_REF] Linke | Pressure-robustness and discrete Helmholtz projectors in mixed finite element methods for the incompressible Navier-Stokes equations[END_REF][START_REF] Linke | Pressure-robustness in quasi-optimal a priori estimates for the Stokes problem[END_REF][START_REF] Lederer | Refined a posteriori error estimation for classical and pressure-robust Stokes finite element methods[END_REF]) if only the pressure (and not the velocity) changes when the external forces acting on the system are modified by a gradient. This property is only valid if there are no harmonic 2-forms, i.e. if the domain is simply connected.

Every vector field f ∈ L 2 (Ω) can be written as ∇ × g + ∇p for some fields g and p. In a bounded domain we only have uniqueness with correct boundary conditions on g and p. As long as these boundary conditions match the ones given in the complex, we have by viewing f as a 2-form, ∇ × g = P B f and ∇p = P B * f . Let f ∈ L 2 (Ω) be another source term and write (ω h , u h , p h , φ h ) (resp. (ω h , ūh , ph , φh )) the solution of (4.14) (resp. (4.14)) for the external force (f, 0) (resp. ( f , 0)).

Theorem 5.3. If f and f are such that ∃p ∈ V 3 , f = f + ∇p then ω h = ωh and u h = ūh . Proof. If there is p ∈ V 3 such that f + ∇p = f then P B (f -f ) = 0 so ∀g ∈ V 1 , f -f , ∇ × g = 0 and in particular ∀g h ∈ V 1 h ⊂ V 1 , f -f , ∇ × g h = 0 .
Therefore P B h (f -f ) = 0 and since we assumed that there were no harmonic 2forms,

P h (f -f ) = P B * h (f -f ). Thus we can find ξ h ∈ V 3 h such that ξ h ⊥ H h and d * h ξ h = -P h (f -f ). Moreover (0, 0, ξ h , 0) verifies: ∀(τ h , v h , q h , χ h ) ∈ V 1 h ×V 2 h ×V 3 h ×H h : 0, τ h -0, ∇ × τ h = 0, νd0 + l 3 0 + l 5 0, v h -ξ h , ∇ • v h ) = f -f , v h , ∇ • 0, q h + 0, q h = 0, ξ h , χ h = 0.
Therefore by linearity and uniqueness of the solution we have (ω h , u h , p h , φ h ) = (ω h , ūh , ph , φh ) + (0, 0, ξ h , 0) and ω h = ωh , u h = ūh . )). However since both ω h and u h remain unchanged when the external force is modified by a gradient, as long as u 0 h = ū0 h and ω 0 h = ω0 h an immediate recursion shows that for all step n, u n h = ūn h and ω n h = ωn h . 5.4. Pressure independent estimate. The pressure-robustness allows us to remove the pressure from the error estimate.

Theorem 5.5 (Pressure-robust estimate). Let (ω, u, p, φ) = (u 1 , u 2 , u 3 , u p ) be the solution of the continuous problem (4.8) and (ω h , u h , p h , φ h ) = (u 1h , u 2h , u 3h , u ph ) be the solution of the discrete problem (4.14). Then it holds:

ω -ω h + ∇ × (ω -ω h ) + u -u h inf χ∈V 1 h ω -χ V 1 + inf v∈V 2 h u -v V 2 .
Proof. Consider an alternative problem with the same l 3 and l 5 as before but with the source term replaced by f2 := f 2 -∇p. Let (ω, ũ, p, φ) and (ω h , ũh , ph , φh ) be respectively the continuous and discrete solution to this alternative problem. By construction, we have p = 0, φ = 0. Hence the estimate of Corollary 4.17 gives: ω -ωh + ∇ × (ω -ωh ) + ũ -ũh inf

χ∈V 1 h ω -χ V 1 + inf v∈V 2 h ũ -v V 2 .
We conclude since by the pressure-robustness we must have: ω = ω, ũ = u, ωh = ω h , ũh = u h .

6. Numerical simulations. We give the results of three numerical simulations done to demonstrate the validity of our scheme. We checked the norm of the divergence angular velocities Ω i and Ω o respectively, and both of height a. The system is closed by two fixed lids at the bottom and top ends. We characterized the system by two geometric parameters: η = R i /R o and Λ = a/d with the gap d = R o -R i . We also need to define two quantities: the inner Reynolds number Re i = Ω i R i d/ν and the outer Reynolds number Re 0 = Ω o R o d/ν where ν is the kinematic viscosity. It is a well known fact that (for an infinite height a) at low speed the flow is steady and fully azimuthal and that vortices start to form at a critical speed. Since a is finite we expect to see vortices near the lids for speeds way under the critical value (they are however fundamentally different from the Taylor vortices, see [START_REF] Ch | Spiral vortices traveling between two rotating defects in the Taylor-Couette system[END_REF]). We compare results obtained from our code with the reference [START_REF] Gebhardt | The Taylor-Couette eigenvalue problem with independently rotating cylinders[END_REF][START_REF] Ch | Spiral vortices traveling between two rotating defects in the Taylor-Couette system[END_REF]. The simulations are done starting from a fluid at rest with no slip boundary condition and studied at t = 0.1. The time step depends on the angular velocities but is most often at δt = 0.001 and only decreased for the highest velocities. We check the value of Re i at which the transition occurs for various values of Re o and η. We see very good agreement, though we used a much coarser mesh and a smaller aspect ratio Λ of 10 instead of 20, for computational cost reasons. We display some values taken on the half plane y = 0, x > 0.

Fig. 2 . 1 .

 21 Fig. 2.1. Degrees of freedom on reference elements for the two lowest polynomial degrees r = 1 and r = 2.

Remark 5 . 4 . 1 h and ω n- 1 h(more precisely l n 3 = l 3 (u n- 1 h) and l n 5 =

 54113315 Consider a time iterating scheme. At any time step n, the linear maps l n 3 , l n 5 involved in the determination of (ω n h , u n h , p n h , φ n h ) are function of u n-l 5 (ω n-1 h

= 2 , d = 1 Fig. 6 . 1 .

 2161 Fig. 6.1. Convergence rate of the velocity with relative H(div) error on a log log scale.

Fig. 6 . 2 .

 62 Fig. 6.2. Isolines for azimuthal velocity and the pressure for Reo = 0 at Re i = 50 (a) and Re i = 72 (b).

Fig. 6 . 3 .

 63 Fig. 6.3. Stream function on cross section y = 0, x > 0 for η = 0.5, Reo = 0 at Re i = 50 (top) and Re i = 72 (bottom).

Fig. 6 . 4 .

 64 Fig. 6.4. Stream function on cross section y = 0, x > 0 for η = 0.8, Reo = 0 at Re i = 90 (top), Re i = 95 and Re i = 100 (bottom).

Figure 6 .

 6 2 shows a comparison between the azimuthal velocity and the pressure for two values of Re i at η = 0.5 and Re o = 0. In Figure 6.3 we see the stream function (the azimuthal component of the vector potential) on the same system. Stream functions for three values of Re i at η = 0.8 are shown in Figure 6.4. Lastly Figures 6.5 and 6.6 give the aforementioned comparison, our results being shown in red and the reference curve in black.

Fig. 6 . 5 .

 65 Fig. 6.5. Comparison of the critical value of Re i at Reo = 0 for various η.

Fig. 6 . 6 .

 66 Fig. 6.6. Comparison of the critical value of Re i at η = 0.5 for various Reo.

Acknowledgements. I thank Pascal Azerad for his comments and helpful advice. I also thank Kaibo Hu for his comments and suggestions and the anonymous referees for many helpful remarks.

References.

Proof. The same proof as the one of Theorem 4.1 works. The only differences will be a sign in the chosen (v, q), l 5 and l 3 substituted by l * 5 and l * 3 and (l 3 d * u, v) changed to (l * 3 u, d * v), which does not add any difficulty in the proof. Remark 4.4. The proof Theorem 4.1 requires taking λ to be sufficiently large, however for any f 2 ∈ V * 2 , f 3 ∈ (W 3 ∩ H 3⊥ ) , λ 0 > 0, λ 1 > 0 we have the following equivalence:

Therefore, if L λ0 is a bounded isomorphism for a given λ 0 , it easily follows that L λ1 is a bounded isomorphism for any λ 1 > 0, in particular for λ 1 = 1. The same argument works for L λ .

From here onward, we assume that L 1 and L 1 are bounded isomorphisms.

4.2.

Well-posedness of the continuous mixed formulation. As we did in the unperturbed case, we introduce an auxiliary variable in the problem. In the following, u is a shortcut for (u 1 , u 2 , u 3 , u ).

We define B by:

The mixed formulation is:

For (u, p) solution of (4.2), it immediately appears that (d * u, u, p, 0) solves (4.8). Now for (u 1 , u 2 , u 3 , u p ) solution of 4.8, the first line implies that

The second line implies that u 3 ∈ V * 3 . And the last implies that u 3 ⊥ H 3 . Therefore (u 2 , u 3 ) ∈ D, and it obviously solves (4.2). The proof follows from the Babuška-Lax-Milgram theorem along with Lemma 4.8 and inf-sup condition (4.13). In the following hidden constants only depend on

by assumption (4.5) and by definition of K we have ∀ω ∈ V :

We set z := (ξ, z 2 , z 3 , u 3H ). Applying (4.9) we have:

Lemma 4.7. For all u ∈ V , there exists z ∈ V such that z V u V and B(u, z)

Proof. Let c = max( l 3 , l 5 ). We begin with some preliminary computations:

(4.10)

Clearly it is possible to construct a suitable z from a linear combination of (4.10), (4.11) and (4.12). Bounds on norms are easily checked, for example:

Combining the two preceding lemmas gives

of the fluid velocity at every time step, and we found a value vanishing up to machine precision. The first simulation aims to verify the pressure-robustness property, the second is based on an exact and fully 3D solution of the Navier-Stokes equation constructed by Ethier [START_REF] Ethier | Exact fully 3D Navier-Stokes solutions for benchmarking[END_REF]. We use it to check the convergence rate in space, first on a steady problem then on an unsteady problem. The last simulation focuses on a system of two rotating cylinders and shows the good agreement with the theory on the value of the critical speed of the inner cylinder at which Taylor vortices appear. It is based upon [START_REF] Ait-Moussa | Numerical Simulations of Coand Counter-Taylor-Couette Flows: Influence of the Cavity Radius Ratio on the Appearance of Taylor Vortices[END_REF][START_REF] Gebhardt | The Taylor-Couette eigenvalue problem with independently rotating cylinders[END_REF][START_REF] Ch | Spiral vortices traveling between two rotating defects in the Taylor-Couette system[END_REF]. In any case, we took a unit kinematic viscosity and polynomials of degree 2. Our codes are written with the FEniCS computing platform, version 2019.1.0 (See fenicsproject.org and [START_REF] Logg | Automated Solution of Differential Equations by the Finite Element Method[END_REF]) and are available at https: //github.com/mlhanot/Navier-Stokes-feec.

6.1. Pressure robustness. We wish to verify that if the external forces acting on two flows differ only by a gradient, then only the pressure differs between the flows. We took the Stokes no-flow problem in a (3-dimensional) glass (see [START_REF] Gauger | On high-order pressure-robust space discretisations, their advantages for incompressible high Reynolds number generalised Beltrami flows and beyond, en[END_REF]). The setup is rather simple, the mesh is a cylinder along the z axis of height 2.0, base radius 1.0 and top radius 1.5 and the force f derives from a potential:

We start from a fluid at rest and enforce a no slip condition on the whole boundary. In every case we found a velocity equal to zero at the order of the machine accuracy. This is not trivial as the same test conducted with Taylor-Hood elements (P 2 /P 1 ) gave a velocity of norm up to 3.3 × 10 -4 for γ = 7.

6.2. Convergence rate to an exact solution. We have conducted a convergence analysis with an exact solution. The expression for the solution is given by Ethier [START_REF] Ethier | Exact fully 3D Navier-Stokes solutions for benchmarking[END_REF] and depends on two real parameters a and b. It is given by: u =   -a(exp(ax) sin(ay + dz) + exp(az) cos(ax + dy)) exp(-d 2 t) -a(exp(ay) sin(az + dx) + exp(ax) cos(ay + dz)) exp(-d 2 t) -a(exp(az) sin(ax + dy) + exp(ay) cos(az + dx)) exp(-d 2 t)   .

We have performed two sets of experiments: the first with a = 2 and d = 0 and the second with a = 2 and d = 1. The domain consists of a cylinder of height 2 and radius 1. In the latter case computations were done for t between 0 and 1 with a time step of 1.0 × 10 -3 . We set the velocity to be equal to the velocity of the exact solution on the boundary (without enforcing any boundary condition on the vorticity). For the stationary case (d = 0) we start from a fluid at rest, otherwise we start from the exact solution at t = 0. We found a rate of convergence in space for the velocity of order 2.0 in both cases which is in agreement with the theory. Figure 6.1 shows the convergence of the velocity in the relative H(div)-norm

(which is the same as the L 2 -norm since the solution is exactly divergence free) with a log-log scale.

Remark 6.1. Despite the fact that the time discretization is Euler implicit and only first order by taking the time step small enough we can observe the second order convergence rate in space. 6.3. Taylor-Couette flow. This test focuses on Taylor-Couette flow. We follow the work of Gebhardt et al. [START_REF] Ait-Moussa | Numerical Simulations of Coand Counter-Taylor-Couette Flows: Influence of the Cavity Radius Ratio on the Appearance of Taylor Vortices[END_REF][START_REF] Gebhardt | The Taylor-Couette eigenvalue problem with independently rotating cylinders[END_REF][START_REF] Ch | Spiral vortices traveling between two rotating defects in the Taylor-Couette system[END_REF]. The geometry consists in two concentric cylinders of constant radius R i for the inner and R o for the outer, rotating at an


Peng Chen

Charles Paillard

Hong Jian Zhao

Jorge Íñiguez

L Bellaiche

Deterministic control of ferroelectric polarization by ultrafast laser pulses

Ultrafast light-matter interactions present a promising route to control ferroelectric polarization at room temperature, which is an exciting idea for designing novel ferroelectric-based devices. One emergent light-induced technique for controlling polarization consists in anharmonically driving a high-frequency phonon mode through its coupling to the polarization. A step towards such control has been recently accomplished, but the polarization has been reported to be only partially reversed and for a short lapse of time. Such transient partial reversal is not currently understood, and it is presently unclear if a full control of polarization, by, e.g., fully reversing it or even making it adopt different directions (thus inducing structural phase transitions), can be achieved by activating the high-frequency phonon mode via terahertz pulse stimuli. Here, by means of realistic simulations of a prototypical ferroelectric, we reveal and explain (1) why a transient partial reversal has been observed, and (2) how to deterministically control the ferroelectric polarization thanks to these stimuli. Such results can provide guidance for realizing original ultrafast optoferroic devices. 2. Rudner, M. S. & Lindner, N. H.

Introduction

Manipulation of properties of quantum materials utilizing (GHz-THz) high-frequency light is a fascinating topic in modern solid-state physics [START_REF] De La Torre | Nonthermal pathways to ultrafast 653 control in quantum materials[END_REF][2][3][4][5][6] . It has resulted in several breakthroughs in the past decade. Examples include stimulating insulator-metal transitions 7,8 , controlling magnetic domains 9-13 , uncovering hidden phases [14][15][16][17] , inducing superconductivity [START_REF] Budden | Evidence for metastable photo-induced 705 superconductivity in K3C60[END_REF][START_REF] Fausti | Light-induced superconductivity in a 708 stripe-ordered cuprate[END_REF][START_REF] Kennes | Transient superconductivity from electronic squeez-711 ing of optically pumped phonons[END_REF][START_REF] Mitrano | Possible light-induced superconductiv-714 ity in K3C60 at high temperature[END_REF] as well as the formation of photon-dressed topological states [START_REF] Mciver | Light-induced anomalous hall effect 717 in graphene[END_REF][START_REF] Sie | An ultrafast symmetry switch in a weyl 719 semimetal[END_REF][START_REF] Wang | Observation of floquet-bloch states on the surface of a 722 topological insulator[END_REF] .

Light-induced switching of ferroelectric polarization [START_REF] Guo | Recent progress in optical control of ferro-724 electric polarization[END_REF][START_REF] Bilyk | Transient polarization reversal using 726 an intense THz pulse in silicon-doped lead germanate. 727 pss[END_REF][START_REF] Yang | Light-induced reversible control 729 of ferroelectric polarization in BiFeO3[END_REF][START_REF] Li | Optical control of polarization in ferroelectric 732 heterostructures[END_REF][START_REF] Mankowsky | Ultrafast reversal of the ferroelectric polarization. 735[END_REF][START_REF] Grishunin | THz electric field-induced second 737 harmonic generation in inorganic ferroelectric[END_REF][START_REF] Morimoto | Terahertz-field-induced large macro-740 scopic polarization and domain-wall dynamics in an or-741 ganic molecular dielectric[END_REF][START_REF] Chen | Ultrafast terahertz-field-driven ionic 743 response in ferroelectric BaTiO3[END_REF][START_REF] Grübel | Ultrafast x-ray diffraction of a ferro-746 electric soft mode driven by broadband terahertz pulses[END_REF][START_REF] Rana | Understanding the nature of ultra-749 fast polarization dynamics of ferroelectric memory in the 750 multiferroic BiFeO3[END_REF] is also among these most important achievements since it can result in novel optoferroic devices, especially ultrafast nonvolatile ferroelectric memories [START_REF] Scott | Ferroelectric memories[END_REF][START_REF] Guo | Non-volatile memory based on the ferroelec-754 tric photovoltaic effect[END_REF] .

In particular, recent works proposed and demonstrated that the ferroelectric polarization can be reversed by exciting a high-frequency infrared-active phonon mode (to be denoted as auxiliary mode or auxiliary high-frequency mode in the following) that is coupled to the soft mode (which is mostly associated with this polarization), with this effect acting through an intermediate anharmonic driving force [START_REF] Mankowsky | Ultrafast reversal of the ferroelectric polarization. 735[END_REF][START_REF] Först | Nonlinear phononics as an ultrafast route 756 to lattice control[END_REF][START_REF] Subedi | Proposal for ultrafast switching of ferro-758 electrics using midinfrared pulses[END_REF] . Such indirect method stimulates coherent phonon modes and has the potential to achieve polarization switching within a few hundred femtosecond, which will be six orders of magnitude faster than from photovoltaic effects [START_REF] Yang | Light-induced reversible control 729 of ferroelectric polarization in BiFeO3[END_REF][START_REF] Li | Optical control of polarization in ferroelectric 732 heterostructures[END_REF][START_REF] Yi | Mechanism of the switchable photovoltaic effect in 762 ferroelectric BiFeO3[END_REF] . This route was first theoretically proposed by A. Subedi and his coworkers [START_REF] Subedi | Proposal for ultrafast switching of ferro-758 electrics using midinfrared pulses[END_REF] and then partially realized by Mankowsky et al. in their second harmonic generation (SHG) experiment [START_REF] Mankowsky | Ultrafast reversal of the ferroelectric polarization. 735[END_REF] . However, some experimental results contradict predictions from the theory. For instance, the model of Ref. [38] predicted a full reversal of the polarization while mea-57 surements "only" reported a transient reversal [START_REF] Bilyk | Transient polarization reversal using 726 an intense THz pulse in silicon-doped lead germanate. 727 pss[END_REF][START_REF] Mankowsky | Ultrafast reversal of the ferroelectric polarization. 735[END_REF] , with 58 the reversed polarization not even reaching its equilib-59 rium value. It is presently unclear why only a partial 60 transient reversal was observed [START_REF] Mankowsky | Ultrafast reversal of the ferroelectric polarization. 735[END_REF] and whether it is, in 61 fact, possible to achieve a full reversal. Possible rea-62 sons for such paucity of knowledge is that the theory of 63

Ref. [38] may have missed important couplings [START_REF] Mertelj | Comment on "ultrafast re-764 versal of the ferroelectric polarization[END_REF] . More-64 over, the fact that the studied model of Ref. [38] is only 65 one-dimensional in nature also implies that some strik-66 ing features may have been overlooked, especially in sys- Ref. [29], but also provides insight into the light-driven 84 effects. It further reveals and explains how a full rever-85 sal (180 • rotation) can indeed happen in some cases; and also predict a variety of light-induced phase transitions, for which polarization rotates by 60 • , 71 • , 90 • , and 109 • , as a result of a mechanism we coin here as "squeezing".

The discovery of such "squeezing" mechanism further allows us to design a strategy for an ultrafast deterministic control of the polarization and even realize full-reversal in a deterministic manner. that is, from cubic to tetragonal (for which a large P z coexists with a small Q z ), then from tetragonal to orthorhombic (where both P y and Q y get activated and become equal to P z and Q z , respectively) and finally from orthorhombic to rhombohedral (for which we have now large P x = P y = P z and small

Results

As

Q x = Q y = Q z ),
when cooling down the system. Note that P and Q are of opposite sign for the tetragonal (T), orthorhombic (O) and rhombohedral (R) phases, as a direct consequence of the positive signs of some coupling parameters indicated in Table I ±3 µC/cm 2 and 30 ±2 µC/cm 2 that one can find at at 147 230 K, 370 K and 473 K in these three phases, respec-148 tively, in Refs. [42, 43]. However as is often the case with The polarization can indeed be fully reversed within a R-

250
phase by a field of 12.46 MV/cm magnitude (red curve).

251

However, not like previously expected [START_REF] Mankowsky | Ultrafast reversal of the ferroelectric polarization. 735[END_REF][START_REF] Subedi | Proposal for ultrafast switching of ferro-758 electrics using midinfrared pulses[END_REF] , the full re-252 versal does not happen during the pulse but after. In

326

To understand such squeezing effect, we need to have a detailed look at Table I and realize that, for the case of a T-phase having a z-component of the polarization, the following couplings are relevant: the bi-quadratic

Q 2 x + P 2 y Q 2 y + P 2 z Q 2 z ) -2.13 × 10 12 Λ22(P 2 x Q 2 y + P 2 x Q 2 z + P 2 y Q 2 z +P 2 y Q 2 x + P 2 z Q 2 x + P 2 z Q 2 y ) 2.
U (z) 22 =Λ 2 P 2 z Q 2 z (1) U (x-z) 22 =Λ 22 (P 2 x + P 2 y )Q 2 z (2)
When a mid-infrared laser pulse polarized along z is ap- If that is the case, we will thus have achieved a fully de-

448

An equilibrium O-phase, having P x = -P y > 0, along perature (300 K) as the initial state. Three pulses are applied along different directions and at different times, 452 along with a bias dc gate field constantly applied along 453 the [ 110] direction. The magnitude of this dc field is form of eq. 2, but for P z and Q y ) will also contribute to 484 the development of P z . Note also that, since no dc field 485 is along the z-axis, the newly-formed P z can also be neg-486 ative, which corresponds to the channel 2 in fig. 3 (a) 487 and which we also observed in our other numerical ex-488 periments (see sec. V of SM). We have thus succeeded to tered on Nb ions. In the KNbO 3 cubic structure, the 559 zone-center P mode is soft (in the sense that the square of its frequency is negative at 0 K) while the Q mode is also located at the zone center but possesses a high and positive frequency. The details of the mode frequency and eigenvector for both P and Q can be found in Sec.I of SM. The local vectors corresponding to the inhomogeneous strains are technically centered on K ions. The homogenous strain is defined with respect to cubic symmetry and has six independent components η i in Voigt notation.

The potential energy U tot of H ef f has four main contributions:

U tot = U F E ({P}, {u}, {η}) + U aux ({Q}, {u}, {η}) + U int ({P}, {Q}) + U elastic ({η}) (1)
a purely elastic one U elastic ; a second one related to the ferroelectric soft mode and its interaction with strains U elastic , exactly like in Ref. [START_REF] Zhong | Competing structural instabilities in cubic perovskites[END_REF] for BaTiO 3 ; a third one related to the high-frequency auxiliary mode and its interaction with strains U aux , which has the same analytical form than the second energy but when replacing P by Q; and a fourth one that gathers the direct interactions between the ferroelectric soft mode and high-frequency auxiliary mode U int , and which are indicated in Table I. The details of all these four energyies can be found in Section II of SM.

Such effective Hamiltonian is generated with respect to cubic (C) symmetry and can be involved in phase transitions among sub-group ferroelectric structures, such as tetragonal (T), orthorhombic (O), and rhombohedral (R) phases. Note that this new effective Hamiltonian resolves symmetry-broken issues expressed in previous works [START_REF] Mankowsky | Ultrafast reversal of the ferroelectric polarization. 735[END_REF][START_REF] Subedi | Proposal for ultrafast switching of ferro-758 electrics using midinfrared pulses[END_REF][START_REF] Mertelj | Comment on "ultrafast re-764 versal of the ferroelectric polarization[END_REF] , and contains for the first time all the symmetry-allowed forms within fourth order. In addition, for the first time too, the model includes the full spatial degree of freedoms of both P and Q, i.e., their Cartesian components along x, y and z axes (which are along the three < 001 > pseudo-cubic directions). The existence of these Cartesian components (P x , P y , P z ) and

(Q x , Q y , Q z )
can, e.g., allow the ferroelectric polarization to rotate (and not "only" reverse its direction) when the Q mode is activated. The model parameters of the H ef f are first determined from 0 K first-principle calculations but then a few of them are adjusted, in order to better agree with measurements but also to allow the calculations to converge (large coupling coefficients can result in simulations going towards infinite energy). This is why, for instance, the coefficient in front of (the third-order-in-P) P 3

x Q x + P 3 y Q y + P 3 z Q z and (the complex) P 2 x P y Q y + P x P 2 y Q x + P 2 x P z Q z + P x P 2 z Q x + P 2 y P z Q z + P y P 2 z Q y have been reduced from its ab-initio value or why the coeffi-cient in the front of Moreover, the kinetic energy K tot of the H ef f contains three parts written with respect to the velocity v p , v q , and v u that correspond to the order parameters, P, Q, and u:

P x P y Q x Q y +P y P z Q y Q z +P z P x Q z Q x
K tot = N i 1 2 M p v 2 p,i + 1 2 M q v 2 q,i + 1 2 M u v 2 u,i (2)
where the effective masses are defined as M p , M q , and 616 M u respectively. Their values are fitted to produce

  67 tems that can adopt polarization along different crys-68 tallographic directions. For instance, is it possible to 69 make the polarization rotate rather than reverse by ac-70 tivating an auxiliary high-frequency mode? Can we also 71 expect novel effects when applying such light-induced in-72 direct method to different structural phases, each having 73 its own direction for the ferroelectric polarization, such 74 as rhombohedral, orthorhombic and tetragonal states in 75 ferroelectric perovskites? 76 The aim of the present study is to answer all these 77 questions by employing an original atomistic scheme 78 that includes the three-dimensional soft mode, the three-79 dimensional auxiliary high-frequency mode and all their 80 relevant couplings. As we will see, this atomistic ap-81 proach not only reproduces a transient partial reversal 82 analogous to the one observed in the experiments of 83

  Figures 1 (a) and 1 (b) report the temperature behavior of the three Cartesian components of the supercell average of the P and Q modes, respectively. These figures indicate that our model can nicely reproduce the sequence of phase transitions observed in experiments,

149effective

  Hamiltonians 44 , some predicted transition tem-150 peratures are lower than the experimental ones 43 , namely 151 640 K versus 708 K for the C-to-T transition, and 345 152 K versus 498 K for the T-to-O transition. On the other 153 hand, the R-to-O transition temperature is well repro-154 duced here: 260 K for the effective Hamiltonian and 263 155 K in measurements. Interestingly, we also numerically 156 found that, if the Q degrees of freedom are turned off, 157 the O-phase vanishes and only two phase transitions oc-158 cur: one from T-to-C at a lower temperature around 500 159 K and a second one from R-to-T around 400 K. Such 160 fact highlights the importance of incorporating both P 161 and Q in order to accurately model properties of KNbO 3 . 162 Note, however, that a previous model without Q did find 163 an O equilibrium phase, but for some narrow tempera-164 ture range ( 50 K) 45 . One can also easily imagine the 165 development of a model having P as the only polar de-166 gree of freedom, but where the effect of Q would be to 167 renormalize the coefficients associated with P.

  168

FIG. 1 .

 1 FIG. 1. Results to be compared with experiments: P (a) and Q (b) modes as a function of temperature; Ratio involving the total polarization at 240K as a function of time, for partial-excitation (d) and full-excitation (e) starting from a R phase. Panel (c) shares the same horizontal axes with Panels (d) and (e) to indicate the laser pulse shape and time evolution. Pink regions are to indicate the full-width-halfmaximum of the Gaussian enveloped laser pulse.

169 170- 1 )

 1691 fig.1 (d)) until such dip becomes negative for 12.47 215

FIG. 2 .

 2 fig.1 (e) is equal to 0.5 and -0.2 for E=8.90 MV/cm and 266

352

  the decreasing of the polarization component that paral-353 lels to the field should be dominant over the increasing 354 of the perpendicular component. Note that the squeez-355 ing effect decides along which direction the polarization 356 should be induced but not the sign of this polarization. 357 For instance, the final result at 2 ps of fig. 2(c) for the 358 axis for which P is parallel to could also have been [110], 359 [1 10], or [ 11 0]. Note too that we have also performed 360 various numerical experiments such, as, e.g., switching 361 in and off some interactions in the simulations, that con-is also realized that the Λ 3 and Λ 1 couplings in Table I decide the direction of the Q mode with respect to that of the P mode, details can be found in sec. IV of SM. It is worthwhile to emphasize that this squeezing effect is a new phenomenon that originates from the highfrequency auxiliary mode. It should not be accounted as the reduction of the polarization by the thermal effect from the light 51-53 because fig. 2 clearly shows that though one component (P z ) decreased to zero, other components (P x and P y ) condense along with it, while heating will typically result in a reduction of the polarization magnitude overall. This squeezing effect also allows us to understand the results about partial (transient) and full (permanent) reversal of polarization in fig. 1 (d) and (e) for some fields within a R-phase. As a matter of fact, pulses with fields oriented along the [11 1] axis and with the resonant 18 THz frequency induce a strong oscillation of the Q mode along such [11 1] axis within the R-phase. During the large oscillation of Q (and resulting large Q 2 x , Q 2 y and Q 2 z ), the P mode starts to decrease its magnitude towards zero according to the squeezing mechanism along all three Cartesian directions. After the exit of the 200 fs full-width-half-maximum pulse, P has just reached and oscillated around zero with net momentum towards the reverse [ 11 1] direction. If it is a full-excitation case (fig. 1 (e), in which all unit cells are subject to the electric field), P can continue to grow along the [ 11 1] direction and the full reversal of the polarization then occurs. If it is a partial-excitation case (fig. 1 (d), in which 10x10x10 cells within the 12x12x12 supercell are allowed to be coupled to the electric field), P can partially grow along [ 11 1] but then quickly reverses its direction and grows back to [11 1], as a result of the interaction with the unexcited part of the material; the transient partial reversal of the polarization thus happens. We also note that Abalmasov 54 proposed that the depolarizing electric field can yield a similar transient partial reversal in thin films. Note that the Supplemental Materials also explain why the polarization can rotate, rather than being reverted, under full-excitation for some fields. Deterministic full-reversal strategy. As can be seen in fig. 1 (e), the ferroelectric polarization full reversal is sometimes replaced by a polarization rotation. Consequently, a single laser pulse, such as the one used in the experiment of Ref. [29], will not reverse the polarization of a three-dimensional ferroelectric in a deterministic manner. (Presumably, the reversal will not be deterministic for a one-dimensional compound like LiNbO 3 either, unless one implements a way to avoid the back-switching to the original state.) However, the understanding of the squeezing effect can enable us to further design a strategy at room temperature to realize a full-reversal of the ferroelectric polarization by laser pulses aided with a tiny dc and constantly applied gate field, which is of high importance for building novel ferroelectric-based devices.

FIG. 3 .

 3 FIG. 3. Deterministic full-reversal strategy: (a) protocol of the three-pulse reversal process; the yellow background boxes are for the initial and excited polarization phases and the pink background boxes are used to indicate the laser pulses that are along x, y, and z-axis; The labels of the phases P is defined as P + Q to indicate the total polarization; The blue arrows are used to indicate the phase transition channels after the laser pulses; Responses of the Q (b) and P modes (c) excited by a sequence of three "200fs-18THz" pulses activated at 10, 20, and 30 ps and having their polarized direction along [100], [010], and [011] directions, respectively. For each of these pulses, the magnitude of E in Ee -2ln2( t τ ) 2 cos(2πωt) is 6.17 MV/cm.

  439

  440 terministic (3-pulse process) reversal of the polarization 441 from the initial to final states. It is worth to mention that 442 the use of several pulses is a common strategy to switch 443 ferroelectric polarization in deterministic manner 55 and 444 has also been suggested by other theoretical work 56 . 445 In figs. 3 (b) and (c), we present one representative of 446 our many numerical examinations, in order to check our 447 proposed strategy and to explain it in even more details.

489

  fig.3 (a), when changing the magnitude of the laser 523

  f ) is developed for such material. It has the follow-551 ing degrees of freedom: vectors related to the ferroelectric 552 soft mode (P), high-frequency auxiliary mode (Q) and 553 inhomogeneous strain (u) in each 5-atom unit cell, as well 554 as, the homogenous strain (η). Both P and Q modes are 555 infrared-active modes and can couple to external electric 556 fields (via an energy involving a dot product with such 557 field). Their associated local vectors in the H ef f are cen-558

608

  has been modified (to give correct A and E mode split-609 ting hierarchy, see supplementary material fig.2). Note 610 also that determining such coupling coefficients from first 611 principles is not a unique and easy procedure and thus 612 large error bars can be assigned to these coefficients. It is 613 thus not surprising that they then require to be adjusted 614 to allow for a better comparison with observations.

  615

640647

  617the TO mode frequencies as calculated from the DFT 618 phonon frequencies of the rhombohedral phase at zero 619 K[START_REF] Wang | Polarization dependence of born effective charge and dielectric constant inKNbO3[END_REF] , and are listed in TableIof the SM. 620 621 Dynamic Simulations We employed Monte Carlo 622 (MC) and Molecular Dynamic (MD) algorithms on a 623 12×12×12 supercell that contains 8640 atoms. More 624 specifically, we used parallel tempering 58,59 (PT) MC 625 and Nosé-Hoover thermal state 60-62 MD simulations 626 implemented in LINVARIANT 63 . The same Nosé mass 627 for P, Q, and inhomogeneous strains of 10000 a.u. 628 is used. The Nosé mass for the homogenous strains 629 is 1 a.u. in all the simulations. Periodic boundary 630 conditions were adapted. For each temperature, 150,000 631 PTMC sweeps are firstly performed, with the first 632 100,000 steps as thermalization and the subsequent 633 50,000 steps to compute the phase diagram; then MD 634 simulations are initialized with the MC outputs and 635 500,000 thermalization steps are performed before the 636 statistical evaluations; the time interval of 0.1fs is used 637 in the MD simulations. The authors declare that the coefficient data of the 641 effective Hamiltonian are available within the paper and 642 its supplementary Material; The data that support the 643 findings of this study are available from the correspond-644 ing author upon reasonable request. The code that generates and solves the effective Hamil-648 tonian is available at https://github.com/PaulChern/

  and further discussed later on. Examples of predicted values for the total polarization for R, O, and T phases are 54 µC/cm 2 at 10 K, 48 µC/cm 2 at 300 K and 42 µC/cm 2 at 400 K. They compare reasonably 145 well with the experimental values of 42 ±4 µC/cm 2 , 32

	146

TABLE I .

 I On-site couplings between ferroelectric soft mode P and high-frequency auxiliary mode Q. The unit of these fourth-order coupling coefficients is N m 6 /C 4 .

	Coefficients (Λ) Coupling Energies
	2.02 × 10 13	Λ2(P 2 x

Acknowledgements 841

rather weak, namely 0.0154 MV/cm which is only 3.3% is first applied at 10 ps with E=6.17 MV/cm and along The authors declare no competing interests.