A priori error analysis for the Shliomis ferrofluid model without magnetic regularization

1.1. The model and its energy law

• The Shliomis model couples incompressible Navier–Stokes flow to magnetization and magnetostatic equations, with particle rotation represented through magnetic torque and corotation.Its nonlinear relaxation is orthogonal to magnetization and contributes nonnegative spin dissipation to the energy law.Source: §1.1 · p. 2

1.2. Analytical difficulties and magnetic regularization

• The energy law gives m and h only L∞(0, T; L2(Ω)) control, while discrete energy stability alone does not ensure convergence or error control for magnetic nonlinearities.Different discrete spaces for magnetization and magnetic field can also prevent an admissible discrete test function m−χh.Source: §1.2 · p. 2, §1.3 · p. 3

1.3. Numerical approximation: results and limitations

• The two-potential approach of Benavides, Nochetto and Shakipov [8] provides error estimates without diffusion but is restricted to 2D, omits rotational coupling, and faces a difficult boundary constraint.Its analysis also relies on projection-stability properties not fully verified for the concrete spaces considered.Source: Error estimates for the coupled… · p. 4

1.4. Contributions of the present work

• The conforming semidiscretization applies in two and three dimensions, retains rotational coupling and nonlinear magnetic relaxation, and incorporates the magnetic boundary condition naturally.It uses an independent magnetization field, a single scalar magnetic potential, and an exactly divergence-free velocity approximation; no magnetic-induction potential is needed.Source: §1.4 · p. 4, A formulation in two and three… · pp. 4–5
• Consistent stabilization controls the discrete Maxwell residual, while incompressibility and Hessian symmetry yield a cancellation that closes the nonlinear error analysis.The cancellation transfers derivatives from the magnetization error and collects remaining derivatives into the residual controlled by stabilization, avoiding separate bounds on those derivatives.Source: Consistent stabilization and error… · p. 5
• The scheme is globally solvable, and a priori error estimates hold for sufficiently regular solutions under approximation and projection-stability assumptions.Concrete 2D finite element spaces satisfy these assumptions; the discussed 3D estimates remain conditional on velocity-projection stability.Source: Consistent stabilization and error… · p. 5
• An averaged H2-conforming projector is defined on L1(Ω) and has optimal Sobolev approximation and W2,p(Ω)-stability properties for 1 ≤ p ≤ ∞.For the specified Hsieh–Clough–Tocher and Walkington families, its p = ∞ stability supplies the uniform stability used for the auxiliary magnetic potential.Source: An averaged H2-conforming projector · p. 6

1.5. Preliminary notation and identities

• For v = ∇ψ with twice-differentiable ψ, the Hessian ∇2ψ is symmetric.This symmetry is stated as a preliminary identity.Source: §1.5 · p. 7
• In d = 2, curl applies to vector arguments and is retained for scalar arguments, with planar products preserving the integration-by-parts and cancellation identities.For example, curl u × m is interpreted as (curl u) × m, and the applied-field condition is curl ha = 0; rotational pairings such as (m×h, curl v) are scalar L2 products.Source: §1.5 · p. 7

• The magnetic potential evolution formulation permits independent testing of magnetization and magnetic field in different finite-dimensional spaces.Source: §2 · p. 7
• The initial potential is the mean-zero solution determined by the magnetostatic and normal-flux relations for the initial magnetization and applied field.Source: §2 · pp. 7–8
• The augmented equation imposes a divergence-free field F with zero normal component on the boundary.Source: §2 · p. 8
• The reformulated system retains the applied field through the boundary condition and includes the velocity, magnetization, and potential boundary conditions.Source: §2 · p. 8
• With the initial compatibility condition, the divergence and normal-flux constraint defects remain zero for all times.Source: §2 · p. 8

• H2-conformity gives globally defined weak Hessians that are symmetric almost everywhere, enabling integration by parts without interior-face jump terms.Source: §3 · p. 9
• Stabilization controls the residual ΔΦ + div M in L2(0,T;L2(Ω)), the interior magnetostatic constraint residual needed for the error analysis.Source: §3 · p. 9
• The discrete energy law includes energy change, dissipation, and nonlinear spin dissipation, and is consistent with the formal energy law.Source: §3 · pp. 10–11
• The pressure is recovered from the momentum residual using the range of the transpose divergence operator, with ambiguity confined to its nullspace.Source: §3 · pp. 11–12
• The scheme has a unique global solution and satisfies an energy bound, while its pressure is unique modulo the nullspace of the transpose of the discrete divergence operator.Source: §3 · pp. 8, 11

4. A priori error analysis

• For a smooth weak solution, the analysis uses linear, bounded, time-independent projections, including a strongly divergence-free velocity projection, to define the discrete error framework.The solution is assumed to have compatible initial data.Source: §4 · p. 12, §4.1 · p. 13

4.1. Discrete error equation

• The discrete error satisfies an energy law with dissipation and a subsequent stability bound, yielding a perturbation estimate for the Galerkin error that controls magnetic spin mismatch.The stability bound uses an integrable time-dependent coefficient.Source: §4.1 · pp. 14, 15, 18
• The nonlinear error analysis exploits a cancellation that relies on exact incompressibility and the symmetry of the Hessian of the magnetic potential.These properties are used when manipulating nonlinear terms in the discrete error stability argument.Source: §4.1 · p. 16

4.2. Bounding consistency errors

• Consistency errors are bounded in terms of projection errors and solution norms, with the resulting constant uniform as β tends to zero.The consistency estimate combines bounds for the momentum, magnetization, and cubic nonlinear terms.Source: §4.2 · pp. 18–19

4.3. Error estimates

• The general a priori error estimate is conditioned on finite-element approximation and projection-stability properties, with constants independent of mesh size h.The spaces and compatible projections are specified later, with one estimate remaining conditional.Source: §4.3 · pp. 19–20
• Theorem 4.12 establishes an a priori estimate for the total discretization error under its stated solution-regularity and approximation assumptions, with a constant independent of h and t.With other parameters fixed, the constant is uniform as β tends to zero if the solution norms and approximation and stability constants remain uniformly bounded.Source: §4.3 · p. 20
• The time-integrated spin-gradient mismatch ∫₀ᵗ ‖m × ∇ϕ − M × ∇Φ‖₂ ds is bounded by Ch^{2k} for every t ∈ [0, T], where k := min{ku, km, kϕ}.Source: §4.3 · pp. 20, 22
• The stated solution regularity makes the initial-data norms finite and suffices for the initial-accuracy condition, so no separate spatial approximation bound for auxiliary projections at t = 0 is needed.Source: §4.3 · p. 22

5.1. Velocity-pressure pair

• The Scott–Vogelius pair uses continuous piecewise velocity polynomials and discontinuous piecewise pressure polynomials; its discrete Leray projection is L2-orthogonal onto the discretely divergence-free space.In two dimensions a Fortin projection provides the stated stability and approximation properties, while a comparable three-dimensional operator remains pending.Source: §5.1 · pp. 22–23

5.2. Magnetization

• Continuous piecewise Lagrange elements for magnetization, paired with the L2 projection and Scott–Zhang interpolant, yield the stated approximation and stability estimates.The Scott–Zhang interpolant is W1,∞-stable and optimal in H1; the L2 projection is optimal in L2.Source: §5.2 · p. 23

5.3. Magnetic potential

• The magnetic-potential construction averages degrees of freedom, including derivative-associated ones, to retain H2-conforming macroelements without requiring pointwise derivatives of the input.The potential space has no essential boundary conditions, so volume averages are used for every degree of freedom and only the zero-mean condition is enforced.Source: §5.3 · p. 23
• For 2 ≤ s ≤ r+1 and j ∈ {0, 1, 2}, the approximation error obeys ∥v − J_hv∥_{j,p} ≤ C h^{s−j}∥v∥_{s,p}.J_h maps L1(Ω) into ˜X, is a projector, and reproduces P_r.Source: §5.3 · p. 24
• The supported potential approximation degree k_φ is 3 for full Walkington, 2 for reduced Walkington and full HCT, and 1 for reduced HCT.These bounds use the stated potential regularity without additional pointwise differentiability assumptions.Source: §5.3 · p. 25

6.1. Convergence test in 2D

• The experiments use manufactured solutions on square and cubic domains, with spatial discretization choices of polynomial degree two for velocity, magnetization, and potential.The 2D test uses β = 10^-5; the 3D test uses β = 0.Source: §6.1 · p. 25, §6.2 · p. 26
• In 2D, experimental rates agree with Theorem 4.12, while the L∞L2 velocity rate is one order higher than guaranteed and the magnetic rates show no such extra order.Whether fourth-order Maxwell stabilization explains the observed quadratic magnetic convergence remains open; the authors report evidence that temporal error does not control those rates.Source: §6.1 · p. 26, Table 1 · p. 26

6.2. Convergence test in 3D

• The 3D test is not fully covered by the theory because velocity does not vanish on the boundary and an appropriate Guzmán–Neilan auxiliary projection with the required W1,∞ stability is unavailable.The nonzero trace introduces a boundary term in the error analysis that the discrete energy law does not control.Source: §6.2 · p. 26
• In 3D, L∞(L2) errors for u, h, and m converge nearly third order, the Maxwell residual about second order, and the velocity L2(H1) rate stays near 1.75.The authors note that the velocity seminorm rate may reflect a preasymptotic regime.Source: Table 2 · p. 27, §6.2 · p. 27

Declaration of Generative AI Use

• The authors report using generative AI for mathematical exploration, proof work, literature searches, numerical code, result interpretation, and manuscript drafting, while retaining responsibility for the content.They say all mathematical content, references, implementations, and numerical results were checked.Source: Declaration of Generative AI Use · p. 27

Appendix A. Proof of Theorem 5.1

• Uniform local stability of J_h is proved by bounding normalized nodal and dual functions, then applying rescaling and polynomial reproduction over the shape-regular mesh family.The proof separates element size from shape.Source: Appendix A · p. 27

A.1. Local representation and proof strategy

• Shape regularity gives the element patches uniformly bounded overlap and comparable element diameters.Each patch consists of elements sharing a vertex with the target element.Source: §A.1 · p. 27
• The proof combines the patchwise estimate with polynomial reproduction and Bramble–Hilbert-style approximation properties.Source: §A.1 · p. 28
• The local expansion of J_h uses nodal basis functions, with local dual functions defined through a duality identity.The expansion is taken element by element in the local nodal basis.Source: §A.1 · p. 28

A.2. Splitting the contributions of element size and shape

• The map F_T rescales and translates the normalized element T* to T, while derivative degrees of freedom scale according to their order.The factor h_T^-μ_a accounts for whether a degree of freedom uses function values, first derivatives, or second derivatives.Source: §A.2 · p. 28
• Shape regularity bounds the inradius of every diameter-one normalized element away from zero, allowing its shape to be encoded by a nondegenerate affine matrix.G_T maps a fixed reference simplex to T*; F_T then restores the element's size and position.Source: §A.2 · p. 28
• The normalized simplices share a fixed reference subdivision into three subtriangles in 2D or four subtetrahedra in 3D.The shape map carries the reference faces and subsimplices to their corresponding normalized-element counterparts.Source: §A.2 · p. 29
• The closure of the normalized shape matrices is compact and contains only invertible matrices, excluding degenerate limiting shapes.Uniform bounds on the shape maps and their inverses support this shape control.Source: §A.2 · p. 29

A.3. Uniform bounds on normalized elements

• The normalized nodal functions are represented by polynomial coefficients on a fixed reference subdivision, with a degree bound independent of the element and mesh size.The polynomial degree bound is distinct from the reproduction degree of the macroelement space.Source: §A.3 · p. 29
• Continuity of piecewise polynomials across internal faces requires matching values and gradients, while additional homogeneous constraints define the selected macroelement space.The extra constraints vary by element family; examples include reduced HCT edge conditions and Walkington barycenter or face conditions.Source: §A.3 · p. 30
• The nodal coefficients solve a shape-dependent linear system encoding continuity constraints and prescribed degree-of-freedom values.Unisolvence on nondegenerate normalized elements and compactness of the shape family yield continuous dependence of the coefficients on shape.Source: §A.3 · pp. 30–31
• For every normalized element, the nodal functions are uniformly bounded in Lp for 1 ≤ p ≤ ∞ and derivative orders 0 through 2.Compactness also gives uniform bounds on inverse mass matrices, which determine the dual functions.Source: §A.3 · pp. 29, 31

A.4. Rescaling and local stability

• Rescaling cancels the powers associated with derivative order between each nodal function and its averaged coefficient.The change of variables uses ∇F_T = h_T I.Source: §A.4 · p. 31
• Hölder's inequality and the dual-function bound control each averaged coefficient for v in Lp(T_a).Source: §A.4 · p. 32
• The local stability estimate follows by inserting the nodal and dual-function bounds into the local expansion and using the bounded number of local degrees of freedom.Each neighboring element T_a lies in the patch ω_T and has diameter comparable to h_T.Source: §A.4 · p. 32

A.5. Approximation and Sobolev stability

• Polynomial reproduction follows by subtracting one polynomial on ωT that accounts for all incident elements in (A.2).For s−1≤r, each local projection reproduces the patch polynomial, so (J_h v)|T=q|T; the constant is independent of h and T under shape-regularity and local quasi-uniformity.Source: §A.5 · p. 32
• The local estimate (A.24) covers derivative orders j=0, 1, and 2 for integer s from 2 through r+1.Its bound is expressed in terms of the W^{s,p}(ωT) seminorm.Source: §A.5 · p. 32
• Summing local estimates over elements and using bounded patch overlap yields global bounds for ∇^j(v−J_hv), including p=∞ by taking a supremum.Combining derivative orders and using h≤1 proves (5.8); taking s=2 and j=2 then gives (5.9).Source: §A.5 · p. 32