Variational principles in quaternionic analysis with applications to the stationary MHD equations

In this paper we aim to combine tools from variational calculus with modern techniques from quaternionic analysis that involve Dirac type operators and related hypercomplex integral operators. The aim is to develop new methods for showing geometry independent explicit global existence and uniqueness criteria as well as new computational methods with special focus to the stationary incompressible viscous magnetohydrodynamic equations. We first show how to specifically apply variational calculus in the quaternionic setting. To this end we explain how the mountain pass theorem can be successfully applied to guarantee the existence of (weak) solutions. To achieve this, the quaternionic integral operator calculus serves as a key ingredient allowing us to apply Schauder's fixed point theorem. The advantage of the approach using Schauder's fixed point theorem is that it is also applicable to large data since it does not require any kind of contraction property. These consideration will allow us to provide explicit iterative algorithms for its numerical solution. Finally to obtain more precise a-priori estimates one can use in the situations dealing with small data the Banach fixed point theorem which then also grants the uniqueness.


Introduction
Nowadays, there are many powerful numerical methods available, for instance well-developed Finite Elements or Boundary Elements Methods, that can be used to solve numerically complicated non-linear systems of partial differential equations arising in modern problems of mathematical physics and engineering in dimension n ≥ 3.According to our knowledge, much less has been developed on the level of analytical methods.To get deeper theoretical inside on the existence, uniqueness, regularity and the structure of solutions analytical methods are very useful, too.
In 1989 K. Gürlebeck and W. Sprößig proposed in their first book [24] a new analytic toolkit for the treatment of three-dimensional non-linear elliptic boundary value problems, including the stationary viscous Navier-Stokes system.The main ingredients are quaternionic integral operators that arise from a higher dimensional generalization of complex function theory in the sense of the Riemann approach which considers null-solutions to the Euclidean Dirac operator, see also [15] or the famous volume [36] edited by J. Ryan which addressed particularly applications to singular boundary value operators, including Calderon Zygmund type operators acting on Lipschitz surfaces.During the following two decades M. Shapiro, V.V. Kravchenko, S. Bernstein, P. Cerejeiras, U. Kähler, J. Ryan, F. Sommen, N. Vieira, see for example [7,9,10,12,24,32,33,22], and many others extended this new analytic machinery to also treat non-linear elliptic, parabolic and hypoelliptic systems.
The quaternionic calculus leads to some new explicit criteria for the existence and the uniqueness of the solutions and also provided some information on the regularity behavior.Based on this quaternionic operator calculus also new numerical algorithms have been worked out in the first decade of this century, see for instance [17,21].In [17,21] it was particularly shown how one can directly derive from the continuous quaternionic operator calculus a discrete quaternionic calculus version, which then even guarantees strong convergence instead of weak convergence only.Furthermore, for some particular domains under some special conditions these methods even allowed us to fully explicit representation formulas for the solutions to the Navier-Stokes equations, the Maxwell, Helmholtz and time independent Klein-Gordon equation, see for example [14].
In the series of more recent papers [11,30,31] it has been shown how to apply the classical quaternionic operator calculus even to the much more complicated stationary viscous magnohydrodynamic equations (MHD) that combine the Navier-Stokes system with the Maxwell-equation.See also the more recent book [26] Section 7.7.6.
Since the works of M. Sermagne and R. Temam in 1983 the study of MHD equations have attracted the interested of a constantly and rapidly growing research community -just to mention a few milestones and some more recent developments we recommend the reader for instance to consult [4,5,8,13,16,18,19,20,27,28,29,34,37,38,39,40,44] among many others.
Also complex quaternions and even octonions have been used for instance in [16,41] in the description of the dynamics of dyonic plasmas in an elegant way.
Most of the literature, including our recent papers [11,30,31], in which we presented explicit estimates involving operator norms and the Reynolds numbers under which we obtain global existence and uniqueness of the solutions, use the Banach fixed point algorithm.
We got an explicit geometry independent formula for the Lipschitz constant, however to have the contraction property we can only apply the presented algorithm to small data.This is a general problem when using the Banach fixed point algorithm, see also [6].One way to overcome this problem to also address large data is to work with Schauder's fixed point theorem instead, such as suggested for the Navier-Stokes system also in [6].To apply the ideas of Schauder's fixed point theorem in the framework of the quaternionic operator calculus however some substantial new tools have to be developed first.
In this paper we now aim for the first time according to our knowledge to combine tools from variational calculus with tools from quaternionic analysis.Indeed, our approach then will allow us to develop new methods for showing geometry independent explicit global existence and uniqueness criteria as well as new computational methods using Schauder's fixed point theorem.Again we concretely focus on the stationary incompressible viscous magnetohydrodynamic equations while the instationary case and other cases will be studied in some of our follow-up papers.
Concretely, in Section 2.3 we show how to specifically apply the variational calculus in the quaternionic setting.To this end we explain in Section 2.4 how the mountain pass theorem can be successfully applied to guarantee the existence of (weak) solutions.To achieve this, the quaternionic integral operator calculus serves again as a key ingredient allowing us later in Section 3 to apply Schauder's fixed point theorem.As already pointed out, the great advantage of the approach using Schauder's fixed point theorem is that it is also applicable to large data since it does not require any kind of contraction property.These considerations will allow us to provide explicit iterative algorithms for its numerical solution which we also present in Section 3, see equations ( 35)- (37).A disadvantage is of course that we only get the existence but in general not the uniqueness of a solution by following this approach.For the sake of giving a complete state of the art we also briefly summarize at the end of the paper in the appendix how the quaternionic operator calculus gives a unique solution when the Banach fixed point is applicable such that the reader is provided with an offrounded self-contained presentation.
2. Toolkits for solving the stationary convective MHD equations 2.1.The stationary convective model.We recall the stationary convective case of the MHD equations.To leave it simple we work in the dimension-less setting and consider the system where u = (u 1 , u 2 , u 3 ) represents the velocity of the flow, p the pressure, B = (B 1 , B 2 , B 3 ) the magnetic field, µ 0 is the magnetic permeability of the vacuum, and R e , R m denote the fluid mechanical and the magnetic Reynolds numbers, respectively.
Equations (3) express the incompressibility of the flow and the non-existence of magnetic monopoles, respectively, while (4) indicate the data at the boundary of the domain.From now on, we shall refer to this system as System 1.
This paper addresses the case where the non-linear convective terms (u grad)u, (B grad)u, and (u grad)B are not negligibly small, but where the flow is still viscous.
We remark that all quantities are three-dimensional.Furthermore • and × denote the standard inner product and vector product, respectively.Furthermore, we assume that Ω is a bounded domain in R 3 with a boundary Γ = ∂Ω representable as a finite union of disjoint closed C 2 surfaces (each having finite surface area).More generally, our considerations can also be applied to the case where Γ is a finite union of Liapunov surfaces or even a strongly Lipschitz surface in the sense of [35,36]  The term x 0 =: Sc(q) is called the real part or sometimes the scalar part of the quaternion q and x = x 1 e 1 + x 2 e 2 + x 3 e 3 is its vectorial part.The latter is also denoted by Vec(q).This consideration allows the embedding of R 3 into H by identifying a purely imaginary element x = x 1 e 1 + x 2 e 2 + x 3 e 3 from H with a vector (x 1 , x 2 , x 3 ) from R 3 .Since the product of two quaternions can be expressed in terms of a symmetric and an anti-symmetric part, we obtain the following relation between the symmetric and the anti-symmetric part, as well as between the standard inner product and the vector product, sometimes also called outer product, in R 3 : The quaternionic conjugation is the involutory automorphism • : H → H defined by q = x 0 + x → q = x 0 − x.Using this notation we can express the Euclidean norm in H by where f 0 , f 1 , f 2 , f 3 : Ω ⊂ R 3 → R are its real-valued coefficient functions.Properties such as continuity, etc., are to be understood component-wisely.
In the particular case of pure vectorial functions, i.e. functions with the mapping behavior u = u 1 e 1 + u 2 e 2 + u 3 e 3 : Ω ⊂ R 3 → H the action of D can be expressed as ( 10) Additionally, we remark that ∆u = −D 2 u.
Using the quaternionic System ( 1)-( 4) can be rewritten in the form which will be called System 2 from now on.Here again, R e and R m denote the mechanical and magnetic Reynolds numbers, respectively.Recall here that Vec(DB) = DB since we always have that Sc(DB) = ∇ • B = 0.
Next we define a right (unitary) quaternionic function module as a function space V (over an open domain Ω ⊂ R 3 whose boundary is for instance a smooth Liapunov surface) together with the algebra morphism (also called right multiplication) R : H → End(V ), given as the point-wise multiplication x → (R(a)f )(x) = f (x)a.When V is a Hilbert space we say we have a quaternionic Hilbert module.
The quaternionic Hilbert module H = L 2 (Ω; H) gives rise to the quaternionic valued functional (15) u, v for u, v ∈ H, and where dx = dx 1 dx 2 dx 3 .Associated to this (quaternionic valued) functional is the scalar functional which corresponds to a real-valued inner product and, hence, gives rise to a norm in the classic sense.
Together with this norm and (quaternionic-valued) inner product the set B(Ω; H) := L 2 (Ω; H) ∩ {f : Ω → H| Df (x) = 0, ∀x ∈ Ω} is a right-quaternionic Hilbert module called the quaternionic Bergman module of left monogenic functions.It is a reproducing kernel Hilbert module in the sense that there is a unique kernel B(x, y) satisfying f (x) = B(x, •), f for all f ∈ B(Ω; H).For details, see for instance [15].Associated to it one can consider the quaternionic Bergman projection P : L 2 (Ω; H) → B(Ω; H) given by P f (x) = B(x, •), f .This projection is an orthogonal projection with respect to •, • .
In what follows, we consider the Sobolev space of L 2 (Ω; H) functions such that the functions and its weak derivatives have finite L 2 −norm, together with the space of test functions Also for simplicity sake, we shall denote from now on the spaces H Following [24,25] and others, one has a Hodge type decomposition (in terms of a direct orthogonal sum) of the form We now use the fact that the differential operator D has a left (but not right) inverse.Following for instance [24] and others for all u ∈ C(Ω; H) we have where q 0 (x) = − x |x| 3 .T Ω is called the (quaternionic) Teodorescu transform.Its boundary analogue is the quaternionic Cauchy transform given by Here, dS(y) is the surface element and n(y) denotes the exterior normal unit at y ∈ ∂Ω = Γ.Again, we may recall from [24] the following important properties For our needs it is important to mention that these two relations remain valid in H 1 (Ω; H), cf.[24].Now apply the operator T Ω QT Ω to System 2. The properties (17), (18), application of the quaternionic Cauchy integral formula QF Γ Du = 0, QF Γ p = 0, QF Γ u = 0 (since u| Γ = 0) and the property that Du ∈ Vec(Q) which implies that T Ω QDu = T Ω Du allow us to rewrite System 2 in the following form, cf.
for details [30,31]: This system will be called System 3 in what follows.Also here and in all that is going to follow the symbols R e and R m stand for the dimension-less mechanical and magnetic Reynolds numbers, respectively.These representation formulas, which are deduced by quaternionic function theory tools will play a crucial role our proof of existence where we want to use Schauder's fixed point theorem in order to also admit the case of large data, for its motivation see also [6].When it is clear which domain Ω is considered, the indices Ω and Γ will be omitted in the sequel for the sake of simplicity.

Preparatory norm estimates.
For our specific purposes we need some particular norm estimates involving some the quaternionic operators that we introduced before.
To start, as an application of the identity DT Ω f = f and the orthogonality u, Dp L 2 we may directly obtain the following generalization of Lemma 4.1 from [10]: (Ω; H) be solutions of system 2. Then we have where M (u) := Sc(uD)u − 1 µ 0 Vec((DB) • B).We also require Lemma 2. Let Ω ⊂ R 3 be a domain with the general requirements mentioned at the end of Section 2.1.Then, relying on [24], we have the following estimate on the operator norm (24) T , where λ min (−∆) denotes the minimal eigenvalue of the negative Laplacian −∆.
Proof.This inequality follows from the norm estimate and from the other Furthermore, applying analogous arguments as in [10,24] we can infer Lemma 3. (Norm estimates) Let u, B be a solution of System 2. Then there is a global constant C S > 0 (called Poisson constant) such that Proof.For the first estimate we refer to [9], Lemma 4.1.The second estimate is analogously to the first and the last is simply the boundedness of the Dirac operator from H 1 (Ω) to L 2 (Ω).
Remark 1.Since the T Ω -operator is a bounded operator from L q (Ω), 1 < q < 3/2, to L 2 (Ω) the above statement means that we have norm estimates like Variational formulation.Now we establish the variational formulation for System 2. The usual variational principles establish In the case of the Dirac operator, we get as e i = −e i .Therefore, the Dirac operator is selfadjoint and has real eigenvalues.
First of all, we observe that so that Equation (11) in System 2 becomes In a similar way, we obtain for Equation (12) the variational formulation Moreover, due to Equation ( 13) one restricts the space of functions further as Hence, the variational form of System 2 is the following: find u, B ∈ H 1,div (Ω; H), and p ∈ L 2 (Ω; H), such that, for all v, w ∈ 2.4.Energy functional.From the above variational formulation we obtain, for v = u and w = B, the energy functional with u, B ∈ H 1,div (Ω; H) and p ∈ L 2 (Ω; H).Notice that we are working in the dimensionless setting; otherwise one has to involve a constant correcting the physical units.Due to Equation ( 13) we obtain since the divergence of u is zero.This implies that Dp is orthogonal to the velocity of the flow u and, furthermore, the energy functional is independent on p.
In a similar way, again due to the divergence of u equal zero.This leads to J(u, B, p) = J(u, B) and hence In the first place notice that J is a lower semi-continuous function.It is also easy to see that J is Fréchet differentiable.This allows us to apply the Mountain Pass Theorem which we recall here for convenience, cf. for example [45]: Lemma 4. (Mountain Pass Theorem).Let X be a Banach space, J : X → R a C 1 -functional satisfying the Palais-Smale condition and J(0) = 0. Suppose a) there exists ρ, α > 0 such that J| ∂Bρ ≥ α where B ρ denotes the ball of radius ρ centered at 0, b) there exists e ∈ X \ B ρ with J(e) ≤ 0.
Then J has a critical value c ≥ α and c is given by where Using this Mountain Pass Theorem we are able to establish the main theorem of this section: Theorem 1.The energy functional J(u, B) has at least one minimum.
Proof.Let us consider the energy functional J(u, B) Now consider (u, B).Again note that p ∈ L 2 (Ω; H) has already been eliminated.The boundary of the ball of radius ρ thus is equivalent to For subsequential purposes we note that we always can rely on ( 30) To be able to apply the mountain pass theorem we have to show that Now, equation ( 29) is equivalent to we next use the identity Therefore we can estimate the expression on the left-hand side of inequality ( 31) by The third term is positive and hence it can be dropped.Therefore, we finally arrived at Using that B H 1 ≤ ρ we obtain that This means that for ρ satisfying ρ < In fact, we here get a condition on how large the radius ρ may be chosen.
Choosing (u 0 , B 0 ) sufficiently far away from the origin we may find J(u 0 , B 0 ) < 0 and the conditions for the mountain pass theorem are fulfilled.This means that the functional J(u, B) has at least one minimum.
We also get the following statement.
where we have C u being the coercivity constant over the domain Ω, i.e.
It is easily seen that G is a continuous mapping.Schauder's fixed point theorem ensures the existence of a fixed point (ũ, B) = (u, B) if we guarantee that taken (ũ, B) from a bounded convex set K we have that G(K) is compact.The map's continuity ensures the boundedness of G(K) while proving (u, B) ∈ H 1+ǫ , ǫ > 0, implies compactness of G(K) due to the compact embedding of H 1+ǫ in H 1 .
To prove that such a bounded convex set K does exist we need the equations in the T QT -form, but the fixed point iteration is constructed as a modification and in a different way from that considered in [31].
For simplicity we here assume that h = 0, otherwise one adds the term F ∂K h + T K P K H to the solution for B, where H ∈ H 2 (K) is an extension of h.
Next, the system above can be rewritten in the following way Equations ( 38) and ( 39) can be rewritten in form of Neumann series (−∆) .Then these two conditions are satisfied if ũ H 1 ≤ min{µ/(Re 2 kC D ), 1/(R 2 m kC D ), } where C D is the operator norm of the Dirac operator.
From these representations we may directly infer that µ 0 T QT Sc(ũD) µ 0 T QT Sc(ũD) Next we can estimate and The latter equation actually implies that there is a constant C > 0 such that Since we have Vec(D B)B L q < +∞, for 1 < q < 3/2, from Theorem 6.1.in [3] we have B ∈ H 1+ǫ , ǫ > 0. This ensures the necessary embedding of (u, B) ∈ H 1+ǫ , for small enough ǫ > 0, and concludes the compactness argument for G(K).
We just have proved our main result of this section: (38) and ( 39) has a solution.We remark that the condition of the theorem also implies the boundedness of B although the theory does not provide uniqueness.To get also uniqueness we have to apply stronger fixed point theorems such as the Banach fixed point theorem which is summarized in the following appendix providing an offrounded presentation.

Appendix: Application of Banach's fixed point algorithm
One advantage of using Schauder's fixed point theorem is that it can be applied to large data.A disadvantage consists in the fact that it usually does not guarantee the uniqueness of a solution.To achieve this one needs to use for instance the Banach fixed point theorem, which however can only be applied to small data as a consequence of the contraction property that is required.In some preceding papers, see [30,31,11], we already showed that the quaternionic operator calculus can also be successfully used in this context.For the sake of completeness we finish this paper by adding a concise summary on the application of Banach's fixed point theorem on the level of these quaternionic integral operators so that the reader gets a rather complete picture of the current state of the art.
To start we wish to emphasize that the representation formulas for the velocity u, the pressure p and the magnetic field B presented in System 3 ( 19)-( 22) also hold when the convective terms are high in comparison with the viscous terms.In the case where the viscous terms dominate the convective terms (the exact conditions will be given Theorem 3 and Theorem 4) we may use the following Banach fixed point algorithm to compute u, B and p, see [31]: , where n = 1, 2, . . .Remark 3. Note that whenever the equations ( 19) and ( 20) have a unique solution, then also (21) will have a unique solution.
In practice we can first compute B n by applying the inner iteration: In [31] we proved Theorem 3. Suppose that u n ∈ • H 1 (G) and that u fulfils the condition has a uniquely solution in H 1 (G) .The sequence defined in (40) converges in H 1 (G) to a unique solution of (42).
Sketch of the proof (for all the details, see [31]).Note that to show the convergence of the inner iteration proposed in (40), we can directly use the estimates from Lemma 2 and Lemma 3. From the simple fact that Sc(a ± b) = Sc(a) ± Sc(b), we may infer from Lemma 3 that Sc(B involving again the same positive constant C S .The norm estimates of Lemma 2 and Lemma 3 yield The classical Banach's fixed point theorem now tells us that B (i) n i = 1, 2, . . .converges under the condition (41) to a unique solution for B n .Finally, since each u n is an element of Now one can further estimate using Lemma 2 and Lemma 3 that For details see again [31].The problem is thus reduced to show the convergence of (u n ) n∈N .
Following further [31]  We can say more.
Theorem 4. (cf.[31]) Under the conditions the proposed fixed point algorithm converges to a unique triple of solutions (u, p, B) of the boundary value problem ( 11)- (15), where p is unique up to a constant when G is bounded, and fully unique when G is unbounded.
For the detailed proof we refer [31].
Remark: The interested reader may also consult [26] Section 7.7.6where the authors also followed basically the same line of arguments using a slightly different notation but established rather analogous estimates based on similar norm constants.The results are qualitatively consistent with those that were obtained previously in [31].
Remarks: This method provides us with an explicit estimate on the Lipschitz contraction constant giving an a priori estimate on how many iterations are needed, cf.[24,25].However, since the Lipschitz constant must be smaller than 1, it can only be applied to small data.A further advantage of the quaternionic operator calculus method is that it always guarantees strong convergence when switching to the discretized versions of these quaternionic operators.It is worthwhile to mention that in many (practical) situations the Teodorescu transform can also be replaced by a simpler primitivation operator [2] whose evaluation requires much less computational steps.This is also a significant input and added value of the quaternionic function theory, since these fnction theoretical concept allow the replacement of a 3D-integral by a one-dimensional primitivation under certain circumstances.Details and extensions to the instationary case will be treated in one of our follow up papers. .

H 1
where• is the usual Sobolev space of the L 2 -function having a first weak Sobolev derivative and with vanishing boundary data.In all that follows we shall denote the complementary part of the Bergman projection by Q := I − P where I is the identity operator.Note that Q one can estimateu n − u n−1 H 1 ≤ 2R 2 e C 1 M (u n−1 ) − M (u n−2 ) H 1 .Furthermore, withC 3 := u n−1 H 1 + u n−2 H 1 and C 4 := B n−1 H 1 + B n−2 H 1 ≤ 2 sup n∈N { B n H 1 } we have M (u n−1 ) − M (u n−2 ) H 1 ≤ C S C 3 u n−1 − u n−2 H 1 F u n−1 − u n−2 H 1This leads to the estimate (43)u n − u n−1 H 1 ≤ L n u n−1 − u n−2 H 1Banach's fixed point theorem tells us that we get convergence if L n < 1.