ERROR ANALYSIS OF AN ADI SPLITTING SCHEME FOR THE INHOMOGENEOUS MAXWELL EQUATIONS

. In this paper we investigate an alternating direction implicit (ADI) time integration scheme for the linear Maxwell equations with currents, charges and conductivity. We show its stability and eﬃciency. The main results establish that the scheme converges in a space similar to H − 1 with order two to the solution of the Maxwell system. Moreover, the divergence conditions in the system are preserved in H − 1 with order one.

1. Introduction. The Maxwell equations are the foundation of the electro magnetic theory and one of the basic PDEs in physics. They form a large coupled system of six time-depending scalar equations in three space dimensions and thus pose considerable difficulties to the numerical treatment already in the linear case. Explicit methods like finite differences on the Yee grid [25] are efficient, but to avoid instabilities one is restricted to small time step sizes, cf. [24]. On the other hand, stable implicit methods for time integration can lead to very large linear systems to be solved in every time step. Around 2000 the very efficient and unconditionally stable alternating direction implicit (ADI) scheme was introduced in [20] and [26] for the Maxwell system on a cuboid with linear isotropic material laws. In this scheme one splits the curl operator into the partial derivatives with a plus and a minus sign, see (1.4), and then applies the implicit-explicit Peaceman-Rachford method to the two subsystems, cf. (1.5). In [20] and [26] it was observed that the resulting implicit steps essentially decouple into one-dimensional problems which makes the algorithm very fast, see also Proposition 4.6 of [13] as well as (4.3) and (4.4) below. There are energy-conserving variants of the ADI splitting, see e.g. [4], [5], [11], [18], not discussed here. We refer to [13], [14] and [15] for further references about the numerical treatment of the Maxwell system. Despite its importance, there exists very little rigorous error analysis of the ADI scheme in the literature, and the available results only cover systems without resistancy, currents and charges. For a variant of the scheme, in [5] error estimates have been shown for solutions in C 6 , see also [4] and [11] for two space dimensions.
The paper [13] (co-authored by one of the present authors) establishes second order convergence in L 2 for the linear Maxwell system on a cuboid with the boundary conditions of a perfect conductor. Here the initial data belong to H 3 and satisfy appropriate compatibility conditions, whereas the coefficients are contained in W 2,3 ∩ W 1,∞ . We stress that the scheme is of second order classically and that the needed degree of regularity in [13] is natural for a splitting in the highest derivatives. As in our paper, the results of [13] are concerned with the time integration on the PDE level and do not treat the space discretization. Based on these and our investigations, we expect that one can develop an error analysis for the full discretization in the future, cf. [14] and [15].
In this work we study the complete Maxwell system with conductivity, currents and charges for Lipschitz coefficients and data in H 2 . Compared to [13], we thus have to modify both the scheme and the functional analytic setting for the Maxwell equations, see (1.4), (1.5) and (2.4), which leads to substantial changes in the analysis. We establish the stability of the scheme in L 2 and H 1 , and that it converges of second order in H −1 , roughly speaking, which is the natural level of regularity for our data. Moreover, the scheme preserves the divergence conditions (1.1c) and (1.2) of the Maxwell system in H −1 up to order τ . To our knowledge, such preservation results have only be shown for C 6 -solutions in the case without charges and in two space dimensions for a related scheme in [4], see also [5] for three space dimensions.
We want to approximate the electric and magnetic fields E(t, x) ∈ R 3 and H(t, x) ∈ R 3 satisfying the linear Maxwell equations in Q, t ≥ 0, (1.1b) div(εE(t)) = ρ(t), div(µH(t)) = 0 in Q, t ≥ 0, (1.1c) E(t) × ν = 0, µH(t) · ν = 0 on ∂Q, t ≥ 0, (1.1d) on the cuboid Q, where ν(x) is the outer unit normal at x ∈ ∂Q. Here the initial fields in (1.1e), the current density J(t, x) ∈ R 3 , the permittivity ε(x) > 0, the permeability µ(x) > 0 and the conductivity σ(x) ≥ 0 are given for x ∈ Q and t ≥ 0. We treat the conditions (1.1d) of a perfectly conducting boundary. As noted in Proposition 2.3, the charge density ρ(t, x) ∈ R depends on the data and (if σ = 0) on the solution via ρ(t) = div(εE(t)) = div(εE 0 ) − t 0 div(σE(s) + J(s)) ds, Throughout, we assume that the material coefficients satisfy For the initial fields and the current density we require regularity of second order and certain compatibility conditions in our theorems. In Section 2 we present the solution theory for (1.1). In presence of conductivity, currents and charges, one has to work in the space X div of fields in L 2 satisfying the magnetic conditions µH · ν = 0 and div(µH) = 0 from (1.1) as well as the regularity div(εE) ∈ L 2 (Q) of the charges, see (2.4). The last condition also enters in the norm of X div . (If σ, J and ρ vanish, it is replaced by the equation div(εE) = 0, see e.g. [13].) The electric boundary condition is included in the domain of the Maxwell operator M from (2.3) governing (1.1a) and (1.1b). It is crucial for our analysis that the domain of the part of M in X div embeds into H 1 (Q) 6 , see Proposition 2.2.
In the ADI method one splits curl = C 1 − C 2 and M = A + B, where we put The domains of A and B are described after (3.1). Let τ > 0, T ≥ 1, and t n = nτ ≤ T for n ∈ N. The (n + 1)-th step of the scheme is given by Here we modify an approach developed in [22] for a different situation. Note that the conductivity σ is included into the maps A and B, whereas the current density J is added to the scheme.
In Section 3 we analyze the operators A and B and their adjoints, showing in particular that they generate contraction semigroups (possibly up to a shift). In L 2 (Q) 6 we can proceed as in [13], but we also have to work in the closed subspace Y of H 1 (Q) 6 equipped with the boundary conditions in (1.1d), which leads to significant new difficulties. Proposition 3.6 then yields the main estimates needed for the error analysis of (1.5). We point out that the domains of the parts A Y and B Y of A and B in Y provide a very convenient framework for the analysis of the ADI scheme, see Sections 4 and 6. In Section 4 we explain the efficiency of the scheme. Based on the results of Section 3, its stability in L 2 and H 1 is shown in Theorem 4.2. These estimates should lead to its unconditional stability independent of the mesh size of a spatial discretization. Moreover, if σ, J and ρ vanish, we obtain a modified energy equality for the scheme.
Our main results are proved in the final two sections. By Theorem 5.1, the error of the scheme is bounded in Y * by cT 2 e κT τ 2 times certain second order norms of E 0 , H 0 and J, where c and κ only depend on the quantities in (1.3). To show this core fact, we adopt arguments from [12] and [22] to derive the (rather lengthy) error formula (5.6). It is formulated in a weak sense with test functions in Y which is needed to work in the present degree of regularity of the data. To estimate the error, one then uses the above mentioned Propositions 2.2 and 3.6. In the final Theorem 6.1 we prove a first order bound in H −1 for the error concerning the discrete analogues of the divergence conditions div(µH) = 0 and (1.2). Besides Proposition 3.6, this proof is based on the surprisingly simple exact formula (6.3) for this error, which nicely follows from the structural properties of the scheme.
Concluding, we discuss possible extensions of our work. For data in H 3 one can establish variants of Theorems 5.1 and 6.1 for the L 2 -norms of the errors. For these results one has to develop a rather intricate regularity theory in H 2 for the Maxwell equation and for the split operators A and B. See our companion paper [9] and the thesis [8], where one can also find numerical experiments. The scheme loses its efficiency for non-isotropic (or nonlinear) material laws so that it does not make sense to study matrix-valued coefficients ε, µ, or σ. On the other hand, numerically one could implement the scheme also on a union of cuboids (e.g., an L-shaped domain). Since such domains are not convex anymore, the domain of the Maxwell operator only embeds into H α (Q) for some α ∈ (1/2, 1), see e.g. [3], which should lead to reduced convergence orders. Further technical difficulties arise since some of the arguments in Section 3 heavily depend on the structure of a cuboid. These questions shall be investigated in a later paper. In [14] and [15], an error analysis was given for the full discretization of the Maxwell system (with σ = 0), using the discontinuous Galerkin method and a locally implicit time integration scheme. We expect that one can treat the full discretization for the ADI scheme combining methods in these and our papers.

2.
Basic results on the Maxwell system. We first collect notation and basic results used throughout this paper. By c we denote a generic constant which may depend only on Q and on the constants from (1.3); i.e., on δ, ε W 1,∞ , µ W 1,∞ , or σ W 1,∞ . We write I for the identity operator and v · w for the Euclidean inner product in R m .
Let X and Y be real Banach spaces. On the intersection X ∩ Y we use the norm z X + z Y . The symbol Y → X means that Y is continuously embedded into X, and X ∼ = Y that they are isomorphic. The duality pairing between X and its dual X * is denoted by x * , x X * ,X for x ∈ X and x * ∈ X * ; and the inner product by (· | ·) X if X is a Hilbert space. In the latter case, a dense embedding Let B(X, Y ) be the space of bounded linear operators from X to Y , and B(X) = B(X, X). The domain D(L) of a linear operator L is always equipped with the graph norm · L of L. If Y → X, then the part L Y of L in Y is given by D(L Y ) = {y ∈ Y ∩ D(L) | Ly ∈ Y } and L Y y = Ly. For two operators L and G in X, the product LG is defined on D(LG) = {x ∈ D(G) | Gx ∈ D(L)}.
Let λ belong to the resolvent set of a closed operator L in X. We occasionally need the extrapolation space X −1 = X L −1 of L; i.e., the completion of X with respect to the norm given by x −1 = (λI −L) −1 x X . One then has a continuous extension L −1 : X → X −1 whose resolvent operators extend those of L. If L generates a C 0 -semigroup T (·) on X, then L −1 is the generator of the semigroup T −1 (·) of extensions to X −1 . This procedure can be iterated, providing L −2 : X −1 → X −2 . If X is reflexive, then X L −1 can be identified with the dual space of D(L * ). See Section V.1.3 in [2] or Section II.5a in [10].
Our analysis of the Maxwell system takes place in the space X = L 2 (Q) 6 with the weighted inner product The square of the induced norm · X is twice the physical energy of the fields (E, H), and because of (1.3) it is equivalent to the usual L 2norm. We further use the Hilbert spaces . Theorems 1 and 2 in Section IX.A.1.2 of [6] provide the following facts. The space of restrictions to Q of test functions on R 3 is dense in H(curl, Q) and H(div, Q). The tangential trace u → (u × ν)| Γ on C(Q) 3 ∩ H 1 (Q) 3 has a unique continuous extension tr t : H(curl, Q) → H −1/2 (Γ) 3 , and H 0 (curl, Q) is the kernel of tr t in H(curl, Q). We also have the integration by parts formula Moreover, Section 2.4 and 2.5 of [21] provide the continuous and surjective trace operator tr : . Its kernel is the space H 1 0 (Q). We also have to deal with cases of partial regularity. For instance, take a function , and thus possesses traces at the rectangles Γ ± 1 = {a ± 1 } × Q 1 whose norms in L 2 (Γ ± 1 ) are bounded by c( f L 2 (Q) + ∂ 1 f L 2 (Q) ). In this way, we obtain trace operators tr Γ ± j and tr Γj for j ∈ {1, 2, 3}. They coincide in L 2 (Γ ± j ), respectively L 2 (Γ j ), with the respective restrictions of tr f if f ∈ H 1 (Q). We usually write u 1 = 0 on Γ 2 instead of tr Γ2 (u 1 ) = 0, and so on. The following lemma will often be used to check boundary conditions.
Proof. We only consider Γ − 1 , j = 1, and k = 2 since the other cases are treated analogously. By the above observations, for a.e. (x 2 , Similarly, . Using the definition of weak derivatives, one checks that , is equal to 1 on [a − 1 + 1/n, a + 1 ], and |χ n | is bounded by cn. We then define f n (x) = χ n (x 1 )f (x) for x ∈ Q. By dominated convergence, the functions ∂ 2 f n = χ n ∂ 2 f converge to ∂ 2 f in L 2 (Q) as n → ∞. In the derivative ∂ 12 f n = χ n ∂ 2 f + χ n ∂ 12 f , the second summand 5690 JOHANNES EILINGHOFF AND ROLAND SCHNAUBELT tends to ∂ 12 f in L 2 (Q) again by Lebesgue's theorem. Let S n be the support of χ n . Using formula (2.2) and Hölder's inequality, we deduce On X we define the Maxwell operator This domain includes the electric boundary condition. To encode the magnetic boundary and divergence conditions in (1.1) and the regularity of the charge density ρ = div(εu), we introduce the subspace The above constraints are understood in H −1 (Q), respectively H −1/2 (Γ). The second equation in (2.4) follows from (1.3) by Remark 3.3 in [13] and because of div(εu) = ∇ε · u + ε div u. In the same way one sees that div v belongs to L 2 (Q) if (u, v) ∈ X div . Equipped with the norm given by 2 L 2 , X div is a Hilbert space as the maps div : L 2 (Q) 3 → H −1 (Q) and tr n : H(div, Q) → H −1/2 (Γ) are continuous.
The part of M in X div is denoted by M div . We actually have for k ∈ N. To show this claim, let (u, v) ∈ D(M ) ∩ X div . As in Proposition 3.5 of [13] (for the case σ = 0) one infers that M (u, v) satisfies the magnetic conditions in X div . Moreover, according to assumption (1.3) the function , and (2.5) is shown for k = 1. By induction, (2.5) follows for all k ∈ N.
The spaces H(curl, Q) and H(div, Q) contain rather irregular L 2 -functions, e.g. from the kernels of curl and div. Nevertheless, their intersection embeds into H 1 (Q) 3 if one assumes that either the tangential or the normal trace is 0. See Theorem 2.17 in [3], for instance. 6 , with a constant only depending on the constants in (1.3). Moreover, in the sense of the traces tr Γj the fields (E, H) ∈ D(M div ) satisfy As noted above, also div E and div H are contained in L 2 (Q). The asserted embedding then follows from Theorem 2.17 of [3]. We thus obtain tr t E = ν × tr E and tr n H = ν · tr H, which yields the second assertion.
We now collect the main properties of the Maxwell operators and solve (1.1). Some of these results are contained in Section XVII.B.4 in [7], for instance. The next proposition is true on any Lipschitz domain Q with the same proof. a) The operators M and M div generate C 0 -semigroups (e tM ) t≥0 on X, respectively (e tM div ) t≥0 on X div . Moreover, e tM div is the restriction of e tM to X div , and we have e tM for t ≥ 0. The charge density in (1.1c) is contained in L 2 (Q) and satisfies Proof. 1) If σ = 0, for instance Proposition 3.5 of [13] shows that M generates a contraction semigroup on X. Using the dissipative perturbation theorem, see Theorem III.2.7 in [10], one can extend this result to the case σ ≥ 0. In the same way one shows that the operator matrix in d) with domain D(M ) is a generator. Because of (2.1) it is a restriction of M * , cf. Proposition 3.5 of [13]. So assertion d) has been shown.
2) We observe that the inhomogeneity ( 1 ε J, 0) satisfies the same assumptions as (J, 0) in part b), respectively c). Let the conditions of b) be true. In view of the magnetic conditions in (1.1c) and (1.1d), we introduce the closed subspace X mag = {(u, v) ∈ X | div(µv) = 0, tr n (µv) = 0} of X. As in Proposition 3.5 of [13], one sees that M maps D(M ) into X mag . Hence, the resolvent (λI − M ) −1 for λ > 0 leaves invariant X mag . The same is true for the operator e tM since it is the strong limit of ( n t ( n t I − A) −1 ) n in X for t > 0, see Corollary III.5.5 of [10]. Due to Duhamel's formula, the solution w then takes values in X mag and thus solves (1.1).
in H −1 (Q). This formula leads to the second part of (2.7), and b) is established.
3) For the remaining assertions in a), we take J = 0. Since e tM is a contraction in X, we have w(t) X ≤ w 0 X . From (2.7) we then deduce the bound div(εE(t)) L 2 ≤ div(εE 0 ) L 2 + ct (E 0 , H 0 ) L 2 and that div(εE(t)) tends to div(εE 0 ) in L 2 (Q) as t → 0. Thus, e tM possesses a restriction e tM div to X div , which forms a C 0 -semigroup there and is bounded by c(1 + t). By Paragraph II.2.3 of [10] it is generated by M div .
4) Under the assumptions of part c), we can differentiate (2.6) in X div twice in t (after the substitution r = t − s.) Hence, w belongs to C 2 ([0, T ], X div ) and 0). Inverting M div − I, we thus obtain w ∈ C 1 ([0, T ], D(M div )). Similarly, the equation w = M div w − 1 ε (J, 0) then implies that w is contained in C([0, T ], D(M 2 div )). These arguments also yield the asserted bound in c).
Remark 2.4. In Proposition 2.3 we also assume that σ = 0 or σ ≥ σ 0 for a constant σ 0 > 0. An inspection of the above proof then shows that we can omit the factors (1 + t) in part a) and (1 + T ) in c). If σ ≥ σ 0 , then the constant c also depends on 1/σ 0 .

The split operators. We now analyze the operators
which are the main ingredients of our splitting algorithm. These operators are endowed with the domains Each domain contains one half of the electric boundary conditions in D(M div ), see Proposition 2.2. These traces exist since they fit to the partial derivatives in C 2 u for A and in C 1 u for B. We note that A and B map into X, Neither the divergence conditions nor the magnetic boundary condition for the magnetic field are included in D(A) or D(B). We further write with D(A 0 ) = D(A) and D(B 0 ) = D(B) for the parts without conductivity. As in Section 4.3 in [13], one shows the following basic integration by parts formula. Let u, ϕ ∈ L 2 (Q) 3 satisfy C 1 ϕ ∈ L 2 (Q) 2 , C 2 u ∈ L 2 (Q) 3 , and (For instance, take (u, ϕ) ∈ D(A) or (ϕ, u) ∈ D(B).) We then have   Proof. Lemma 4.3 of [13] says that A 0 and B 0 are skew-adjoint on X, and hence generate a contraction semigroup. This property is inherited by A and B due to (3.2) and the dissipative perturbation Theorem III.2.7 in [10]. In the same way one shows the generator property of the operator matrices defined in part a) on the domains D(A) and D(B), respectively. As in Lemma 4.3 of [13], equation (3.3) implies that these operator matrices are restrictions of A * and B * . They are thus equal to these operators, respectively, since the right half plane belongs to all resolvent sets. The other assertions in b) then easily follow, using also Proposition 2.2.

It follows
For our error analysis we need the restrictions of the above operators to the subspace of H 1 given by We use on Y the weighted inner product with the induced norm · Y . Due to (1.3), this norm is equivalent to the usual one on H 1 . The continuity of the traces implies that Y is a closed subspace of H 1 (Q) 6 . We very often use that maps like (u, v) → (εu, µv) leave invariant Y because of (1.3). Our definitions yield the embedding (3.5) We denote by A Y , B Y , (A * ) Y , and (B * ) Y the parts of A, B, A * , and B * in Y , respectively. Their domains are described in the next lemma. Then Here the respective first line in a) follows from Proposition 3.1 and (1.3). Part b) is a consequence of Lemma 2.1 and it yields the rest of assertion a). In a series of further lemmas we collect the basic properties of the above operators. Let (u, v) ∈ Y . We approximate u 1 =: f and v 1 =: g in Y by maps f n and g n which are the first and fourth components of vectors (u n , v n ) ∈ D(A Y ), respectively. We use smooth cut-off functions χ and are equal to one on [a − j + 1 n , a + j − 1 n ] for n > (2d Q ) −1 =: and j ∈ {1, 2, 3}. Also, ρ (j) n is a standard C ∞ -mollifier with support in [− 1 2n , 1 2n ] that acts on x j . 2) We are looking for functions f n ∈ H 1 (Q) with ∂ 2 f n ∈ H 1 (Q) that converge to f in H 1 (Q) and have zero traces on Γ 2 and Γ 3 . We first set ϕ n = χ n ≥ n 0 . These maps belong to H 1 (Q), vanish near Γ 2 and have trace 0 on Γ 3 . Due to dominated convergence, the functions ϕ n tend to f and χ (2) n ∂ j f to ∂ j f in L 2 (Q) as n → ∞ and for j ∈ {1, 2, 3}. As in the proof of Lemma 2.1, one shows that ((χ (2) n ) f ) converges to 0 in L 2 (Q) since f vanishes on Γ 2 . Summing up, (ϕ n ) has the limit f in H 1 (Q). We then extend ϕ n by 0 to R 3 and define f n = ρ (2) n * ϕ n | Q for all n > . This function and ∂ 2 f n belong to H 1 (Q), and f n vanishes near Γ 2 . For a map h ∈ H 1 (Q) ∩ C(Q) with h = 0 on Γ 3 it is clear that ρ (2) n * h is also equal to 0 on Γ 3 . By approximation, we thus obtain tr Γ3 f n = 0.
3) We next have to construct functions g n ∈ H 1 (Q) with ∂ 3 g n ∈ H 1 (Q) such that tr Γ1 g n = 0, tr Γ3 ∂ 3 g n = 0 and (g n ) has the limit g in H 1 (Q). Let Φ be the linear and bounded Stein extension operator that maps functions in H k (Q) to functions in H k (R 3 ), see Theorem 5.24 in [1]. Extending g by 0 to R 3 , we set ψ n,m = ρ (2) n * ρ (3) n * Φ(ρ (1) m * (χ (1) m g))| Q Q for all n, m > . These smooth maps and their derivatives vanish near Γ 1 . Let η > 0. Arguing as in step 2), we fix an index m = m(η) > such that By the properties of mollifiers, there also exists a number n = n(η) > with Setting g = ψ n, m , we obtain the inequality In view of the needed boundary condition of ∂ 3 g n , we define for x = (x , x 3 ) ∈ Q and n > . The functions g n and ∂ 3 g n = χ (3) n ∂ 3 g are contained in H 1 (Q). Moreover, the traces of g n on Γ 1 and of ∂ 3 g n on Γ 3 are zero by construction. The integrands of r n and their derivatives with respect to x 1 and x 2 are uniformly bounded by a constant, and these maps tend to 0 pointwise a.e. as n → ∞. By dominated convergence, the functions r n , ∂ 1 r n and ∂ 2 r n thus converge to 0 pointwise a.e. and then in L 2 (Q) as n → ∞. The same is true for n − 1)∂ 3 g. As a result, (g n ) has the limit g in H 1 (Q). The other components of (u, v) are treated in the same way.
We set Proof. Again we only consider A Y . Let (u, v) ∈ D(A Y ). In view of the boundary conditions in Lemma 3.2, integration by parts yields
2) Let j = 2. We are looking for a function w ∈ D(D 2 ) solving (3.8), where we put h := h 1 . To this aim, we abbreviate µ ∂ 2 w for w ∈ D(D 2 ). As in the proof of Lemma 4.3 in [13], where σ = 0 and κ Y = 0, by means of the Lax-Milgram lemma we obtain a unique map w in D(D 2 ) with Lw = h, and w thus satisfies (3.8). Moreover, Theorem VI.2.7 in [16] shows that L is given by the symmetric, closed, positive definite and densely defined bilinear form , and that L is invertible and self-adjoint in X.
3) To prove that w = 0 on Γ 3 , we approximate h by functions h n = χ n h with n > (2d Q ) −1 , cf. step 2) of the proof of Lemma 3.3. As in that proof one sees that h n vanishes if |x 3 − a ± 3 | ≤ 1/(2n) and that (h n ) tends to h in H 1 (Q) as n → ∞. The map w n := L −1 h n ∈ D(D 2 ) thus converges to w in D(∂ 2 ). We take functions φ n ∈ C ∞ c ((a − 3 , a + 3 )) which are equal to 1 on [a − 3 + 1/(2n), a + 3 − 1/(2n)]; i.e., h n = φ n h n . Since L is injective and Lw n = h n = φ n h n = φ n Lw n = L(φ n w n ), we obtain w n = φ n w n = 0 on Γ 3 . Formula (3.9) and assumption (1.3) next yield for k ∈ {1, 2, 3}. The right-hand side tends to 0 as n → ∞, so that w n converges to w in H 1 (Q). As a result, w has the trace 0 on Γ 3 . We set u 1 = w. 4) In a final step we define compare (3.7). The trace v 3 on Γ 3 vanishes since (f, g) ∈ D(A Y ) and ∂ 2 u 1 = 0 on Γ 3 by Lemma 2.1. We differentiate the above equation w.r.t. x 2 and insert the equations Lu 1 = h 1 and (3.8). It follows in L 2 (Q); i.e., u 1 and v 3 satisfy (3.7). This equation also yields that ∂ 2 v 3 belongs to H 1 (Q) and has trace 0 on Γ 2 , as required in D(A Y ). The other components are treated in the same way.
We now easily derive the basic properties of our split operators on Y by means of the Lumer-Phillips theorem, see e.g. Section II.3.b in [10]. Recall the definition of κ Y in (3.6).
Proof. The generation property and the resolvent bounds follow from Lemmas 3.3-3.5 and the Lumer-Phillips theorem. Let 0 < τ < 1 κ Y and h ∈ Y . The function w = (I − τ L Y ) −1 h then satisfies h = (I − τ L Y )w = (I − τ L)w and hence w = (I − τ L) −1 h, which is the asserted restiction property. For Re z > 0 we define It is easily seen that for all 0 < τ ≤ τ and a constant τ ∈ (0, 1/κ Y ) only depending on κ Y . Since L Y − κ Y I generates a contraction semigroup on a Hilbert space, Theorem 11.5 of [17] provides an H ∞ -functional calculus for κ Y I − L Y and the desired estimate for all τ ∈ (0, τ ), where we recall (3.4). (In this argument we use the complexification of Y and the natural extension of L Y .) We set τ 0 = min{ τ , (2κ Y ) −1 }.

4.
The ADI splitting scheme. Let τ > 0. We set t n = nτ for n ∈ N 0 and assume that (J(t), 0) ∈ D(A) for all t ≥ 0. The ADI splitting scheme is given by the operators For n ∈ N 0 and w 0 ∈ D(B), we further write Proposition 3.6 then yields w n ∈ D(B Y ) and w n+1/2 ∈ D(A Y ) for all n ∈ N 0 .
The operators S (k) τ contain implicit steps. For σ = 0 and J = 0 it was pointed out in [20] and [26] that these steps decouple into (essentially) one dimensional problems, see also [13]. We now extend this observation to our setting. To this aim, for λ ∈ {ε, µ} we define the operators D (1) on the domains D(∂ 22 )×D(∂ 33 )×D(∂ 11 ) and D(∂ 33 )×D(∂ 11 )×D(∂ 22 ), respectively, where D(∂ kk ) is the set of f ∈ L 2 (Q) with ∂ kk f, ∂ k f ∈ L 2 (Q) and f = 0 on Γ k . Let w 0 ∈ D(B Y ) and ( 1 ε J(t), 0) ∈ D(A Y ) for all t ∈ R. Remark 4.1 then yields (E n , H n ) ∈ D(B Y ) for each n ∈ N. The above definitions lead to Because of this regularity, we can insert the second equation into the first one and infer in L 2 (Q) 3 , using curl = C 1 − C 2 . Observe that E n+1/2 belongs to the domain of D (1) µ . Due to (4.2), the implicit part of these equations splits into (essentially) one dimensional problems; one only has to solve parameter-dependent elliptic equations in one space variable. Similarly, the second half step is rewritten as 2ε (J(t n ) + J(t n+1 )). As above we then obtain the essentially one-dimensional problem Using the Cayley transforms γ τ (L) from (3.4) and induction, we further deduce from (4.1) the closed expression of the scheme The constants c > 0 only depend on the constants from (1.3).
Remark 4.3. a) In the above result, we can drop the factor 1 ε in the assumptions if ε also belongs to W 2,3 (Q) since H 1 (Q) → L 6 (Q). b) In Theorem 4.2, for σ = 0 and J = 0 the inequality in X is actually an equality since then the operators A and B are skew-adjoint in X by Lemma 4.3 of [13], and hence their Cayley transforms are unitary in X. This can be viewed as a modified energy preservation of the scheme in the conservative case.

5.
Convergence of the ADI scheme. For the semigroups from Proposition 2.3, τ ∈ (0, 1] and j ∈ N 0 , we define Λ j+1 (τ ) = 1 j!τ j+1 τ 0 s j e (τ −s)M ds and Λ div j+1 (τ ) = 1 j!τ j+1 τ 0 s j e (τ −s)M div ds, as well as Λ 0 (τ ) = e τ M and Λ div 0 (τ ) = e τ M div . By Proposition 2.3, these operators are uniformly bounded on X, respectively X div , and Λ div j (τ ) is the restriction of Λ j (τ ) to X div . Standard semigroup theory shows that these operators leave invariant D(M k ), respectively D(M k div ), and commute with M k , respectively M k div . For j ≥ k they actually map into D(M k ), respectively D(M k div ). One can further check that Our first main result establishes the second order convergence of the ADI scheme in Y * . According to (3.6) the number κ Y ≥ 0 only depends on the constants in (1.3), and we have κ Y = 0 in the case of constant coefficients. We use the number τ 0 > 0 from Proposition 3.6, which only depends on κ Y . H) be the solution of (1.1) and w n = S J τ,n · · · S J τ,1 w 0 be its approximation from (4.1). For all nτ ≤ T and (ϕ, ψ) ∈ Y , we then have We insert (5.2) for s = τ into the definition of S J τ,n+1 from (4.1) withw = S J τ,n · · · S J τ,1 w 0 ∈ D(B) and obtain with the remainder In the next step we take the inner product in X of the fields y = (ϕ, ψ) ∈ Y with the difference of (5.4) and (5.3). In the following we write A * Y and B * Y instead of (A * ) Y and (B * ) Y , respectively. Putting several operators as adjoints on the side of y, we arrive at the formula S J τ,n+1 · · · S J τ,1 w 0 − w((n + 1)τ ) (ϕ, ψ) X (5.5) The terms r k (τ ) and R k (τ ) are bounded in X by cτ 2 (k+1)τ kτ (J (s), 0) X ds. Propositions 2.2, 2.3, 3.1 and 3.6 then imply where we use τ ≤ 1 and c only depends on the constants from (1.3). These arguments also yield the addendum. 6. Almost preservation of the divergence conditions in H −1 . The solution (E, H) of (1.1) fulfills the Gaussian laws (2.7) and div(µH(t)) = 0. We now show that the scheme (4.1) satisfies a discrete analogue of these divergence conditions up to an error of order τ in H −1 (Q). We recall that the numbers κ Y ≥ 0 and τ 0 > 0 from (3.6) and Proposition 3.6 only depend on the constants in (1.3), and that κ Y = 0 if the coefficients are constant. In the proof, see (6.4), we will see that the integral on the left-hand side of the above inequality can be replaced by the sum n−1 k=0 τ 2 div(J(t k ) + J(t k+1 )), 0 .
Moreover, as in Remark 4.3 one can drop the factor 1 ε in the assumption if ε also belongs to W 2,3 (Q).
Proof. We first derive a recursion formula for the divergence of w n which is then estimated by means of Propositions 3.1 and 3.6. Remark 4.1 says that w n belongs to D(B Y ) and w n+1/2 to D(A Y ). We often use this regularity, these propositions and that τ ≤ 1 without further notice. Take k ∈ N 0 with k + 1 ≤ T /τ . 1) As before (4.3), the definition (4.1) yields To obtain a convenient recursion formula, we insert the second line into the first one and the first line into the second one. It follows in L 2 (Q) 6 . Using curl = C 1 − C 2 and the definition (4.2), we reorder the above identity and infer the first half of the recursion step in L 2 (Q) 6 . Similarly, (4.1) leads to the expression 2ε (J(t k ) + J(t k+1 )) in Y . Proceeding as in the first half step, we conclude in L 2 (Q) 6 . Again with (4.2) and curl = C 1 − C 2 , this equation implies the second step of the recursion (1 − στ 4ε )(J(t k ) + J(t k+1 )) 0 + τ 2 0 curl τ 2ε (J(t k ) + J(t k+1 ))