Stochastic Modeling and Regularity of the Nonlinear Elliptic Curl-Curl Equation

This paper addresses the nonlinear elliptic curl-curl equation with uncertainties in the material law. It is frequently employed in the numerical evaluation of magnetostatic fields, where the uncertainty is ascribed to the so-called B-H curve. A truncated Karhunen-Lo\`eve approximation of the stochastic B-H curve is presented and analyzed with regard to monotonicity constraints. A stochastic non-linear curl-curl formulation is introduced and numerically approximated by a finite element and collocation method in the deterministic and stochastic variable, respectively. The stochastic regularity is analyzed by a higher order sensitivity analysis. It is shown that, unlike to linear and several nonlinear elliptic problems, the solution is not analytic with respect to the random variables and an algebraic decay of the stochastic error is obtained. Numerical results for both the Karhunen-Lo\`eve expansion and the stochastic curl-curl equation are given for illustration.

1. Introduction. Today, it is increasingly acknowledged that uncertainty quantification is an important part within simulation based design. It allows to control the risk of failure as technical devices are designed and operated closer to their physical limits. When the underlying physics are modelled by partial differential equations, uncertain inputs are typically identified with the material's constitutive relation, geometries, initial data or boundary values. In particular the case of random material coefficients for linear equations has received considerable attention in recent years, see [26,52,9,6] among others. In electromagnetics the coefficients are frequently modeled to be piecewise constant on subdomains. In a stochastic setting, which is adapted here, neglecting anisotropy, on each subdomain the material coefficient can be represented by a single random variable. An important exception is the magnetic properties of ferromagnetic materials, that, in the anhysteretic case, are expressed through a map where H and B are the magnetic field and flux density, respectively. Magnetic saturation effects, incorporated through the nonlinear dependency on the field magnitude, cannot be neglected in many situations and have been found to be sensitive to uncertainties. In this setting the input randomness is modelled by an infinite-dimensional random field and its discretization must be accomplished. This is complicated by the fact, that each trajectory of the material law has to fulfill smoothness and shape requirements. An example of a material law is given in Figure 1 on the left, showing clearly a monotonic behavior of f . An important aspect of this work is the shapeaware modelling of uncertainties in view of a constraint for the derivative (1.2) 0 < α ≤ ∂f (·, s) ∂s ≤ β < ∞, s ≥ 0.
To this end, in the literature, closed-form parametric material models have been employed. Among others, we mention different forms of the Brillouin model, see, e.g., [48,45] or the Brauer model [12]. These models are appealing due to their simplicity and a physical interpretation of the parameters can often be given. Finite dimensional random fields are readily obtained by using random-instead of deterministic parameters. However, there is a lack of flexibility due to the a priori fixed dimensionality and specific shape of the analytical functions used. Moreover, the model parameters have been found to be correlated [45] and the model might not be used directly in stochastic simulations. The (linear) truncated Karhunen-Loève expansion [35,29] is known to be a flexible and efficient tool to approximate random fields with high accuracy and to separate stochastic and deterministic variables. It is applied in this paper in view of a stochastic material law with regularity and shape constraints (1.2). The regularity can be controlled by the smoothness of the covariance function whereas the shape constraints imply restrictions on the truncation order, or alternatively, on the uncertainty magnitude. It is observed that this magnitude is dependent on the correlation length of the process, in accordance with [8, pp. 1281-1283] in the context of uniform coercivity constraints. Numerical examples with data supported from measurements given in [45] support the findings.
Given random material input data, we discuss a stochastic nonlinear magnetostatic formulation. For related work see [48,12] and [18] for the full set of (linear) Maxwell's equations in a stochastic setting. Following a frequently used procedure in the literature [8,39], we first analyze the modelling error arising in magnetic fields through the truncation of the Karhunen-Loève expansion. Then the problem is reformulated as a high dimensional deterministic one and an approximation scheme is presented. The scheme involves linearization, as well as a finite element and collocation approximation in the deterministic and stochastic variable, respectively. We analyze the stochastic regularity to establish the convergence rate of the numerical procedure. A complication arises here due to the specific type of the nonlinearity. In particular the implicit function theorem cannot be applied and no analytic dependency of the solution with respect to the random variables is obtained. Instead, finite differentiability is established using the chain rule and the deterministic regularity of the solution. Numerical examples will complement the findings. Let us also mention, that the problem considered here is related to many other physical problems, e.g., nonlinear heat conduction. In two dimensions the equations reduce to a nonlinear version of the Poisson equation.
The paper is structured as follows. After a brief description of the magnetostatic model problem in Section 2, we introduce randomness in the material law, present the truncated Karhunen-Loève expansion and analyze the respective truncation error in the stochastic problem in Section 3. In Section 4 we will outline the stochastic collocation method and establish its convergence. This involves a higher order sensitivity analysis w.r.t. the stochastic variables. Finally, in Section 5 convergence results will be illustrated by numerical examples. 1.1. Notation. Boldface type is used for vectors u = (u 1 , u 2 , u 3 ) and vectorfunctions. Important function spaces are 3 and u × n = 0, on ∂D}, and for s > 0 see [37]. We write (·, ·) 2 for the L 2 (D) 3 -inner product and · 2 for the associated norm.
The Euclidean norm is denoted by | · |. Also, we introduce the notation g (i) := ∂ i x g (x) for the i-th derivative of a function. For I ⊂ R n , open and bounded, C k Ī denotes the space of k-times differentiable functions with bounded and uniformly continuous derivatives up to order k, endowed with the norm The sets of real non-negative and positive numbers are denoted R + 0 and R + , whereas N 0 , N refer to the sets of natural numbers with and without zero, respectively.
2. The Model Problem. We consider the magnetostatic problem on a domain D, with outer unit normal n and divergence free electric current density div J = 0. We introduce the magnetic vector potential A, such that curl A = B. Then using the material law (1.1), equations (2.1) are transformed into a second order curl-curl problem In (2.2a), ν refers to the magnetic reluctivity, defined by Although ν is the coefficient appearing in the differential equations measurement results are typically given directly for f and thus both will be included in the discussion. Note that, in contrast to f , ν is not necessarily monotonic. From now on we simplify the situation and neglect the spatial dependency in f , i.e., f (x, ·) = f (·), by abuse of notation. An adaption to the important case of piecewise constant material properties is achieved by minor modifications, see also Remark 1. In absence of hysteresis and anisotropy, the nonlinear magnetic material law can be described, following [47,44], by a bijective function Following [44], properties of the deterministic material law f are summarized in the following assumption.
Assumption 1. It holds for the deterministic material law that f is continuously differentiable, (2.5a) As f might have increased differentiability properties, we refer to (2.5a) as minimal regularity assumption. Note, that β can be identified with the reluctivity of vacuum ν 0 . Depending on the problem formulation it might be more convenient to work with the inverse law f −1 . However, in this case similar assumptions can be made.
2.1. The Nonlinear curl-curl Formulation. We proceed with the derivation of a weak formulation of the model problem and a result on existence and uniqueness of a solution. Throughout the paper, we consider a bounded, simply connected polyhedral Lipschitz domain D. A weak formulation of (2.2) relies on (2.7)V = {u ∈ V | (u, grad ϕ) 2 = 0, ∀ϕ ∈ H 1 0 (D)}, the space of functions in V with weak zero divergence. We recall from [32,Corollary 4.4] that the Poincaré-Friedrichs-type inequality holds, for all u ∈V andV can be endowed with the norm u V := curl u 2 , see also [3] for the more general case of multiply connected D. Then for J ∈ L 2 (D) 3 the weak formulation reads, find A ∈V , such that Equation (2.9) can be written more compactly, by introducing the vector function h : Also, letV * denote the dual space ofV , by introducing the operator K :V →V * as we obtain A ∈V as the solution of (2.12) KA, v = (J, v) 2 , ∀v ∈V .
Existence and uniqueness is guaranteed by the Zarantonello Lemma [54] as (2.6) implies see [34], i.e., the strong monotonicity and Lipschitz continuity of the nonlinear operator K. Moreover, we have the estimate (2.14) Remark 1. Typically, for the accurate modelling of magnetic devices an interface problem with several materials, e.g., iron and air, has to be studied. Then, the piecewise defined magnetic reluctivity also satisfies (2.6) and the problem is still found to be well-posed, see [10,31]. Also the extension to multiply connected domains would mainly require a modification of the divergence free condition to ensure the norm equivalence [10]. In the more general case of f (x, ·) with arbitrary x-dependence, f (·, s) must additionally be measurable [53] for all s and consequently 3 + 1-dimensional random fields would occur. However, as our focus lies on the nonlinearity in the stochastic setting, for simplicity, we restrict ourselves to simply connected domains with only one homogeneous nonlinear material.

Uncertainties in the Nonlinear Material Law and Stochastic
Formulation. Randomness is incorporated, as usual, by introducing a probability space (Ω, F, P) and modeling the material law f as a random field f : Ω × R + 0 → R + 0 . A stochastic formulation is based on the following assumption.
Setting for r ∈ R 3 , h (ω, r) = ν (ω, |r|) r, the stochastic curl-curl problem reads a.s. as By Assumption 2 we have a unique solution A ∈ L p (Ω,V ) for all p ∈ N by means of (2.14).

Random Input Discretization by the Truncated Karhunen-Loève
Expansion. In the following we restrict ourselves to an open interval I ⊂ R + and define for f (resp. ν),f := f |Ī . Restricting the uncertainty of the material law to a specific interval is a reasonable assumption in practice and in particular suitable to satisfy the constraints (2.5c) and (2.5d) in the presence of randomness. A globally defined f can be obtained by a differentiable prolongation from I to R + 0 . To be used in computer simulations the input random fieldf : Ω ×Ī → R + has to be discretized. This is achieved here by introducing the (linear) truncated Karhunen-Loève expansion A stronger dependency off M w.r.t. Y n might be obtained by performing a Karhunen-Loève expansion for log(f ), rather than f [6]. However, a linear expansion as (3.2) with a moderate number M , might be beneficial in numerical approximations.
We recall the definition of the expected value and covariance function for s, t ∈Ī. Based on the following assumption, some important properties of (3.2) are briefly recalled here, see, e.g., [35,26,8,49].
Assumption 3. The image of Y n , n = 1, . . . , M is uniformly bounded. Additionally, the covariance satisfies Covf ∈ C l I × I , with l > 2. Consider the self-adjoint and compact operator Tf : L 2 (I) → L 2 (I), given by Then (λ n , b n ) ∞ n=1 are eigenpairs of (3. 6) Tf u = λu, where the (b n ) ∞ n=1 are orthonormal in L 2 (I) and λ 1 ≥ λ 2 ≥ · · · ≥ 0. In (3.2), the random variables (Y n ) ∞ n=1 subject to for λ n > 0, are centered and uncorrelated with unit variance. We assume that they are independent and denote their image with Γ n and set Γ := M n=1 Γ n . Approximation properties of (3.2) are well studied, in particular the L 2 -error is optimal among all M -term approximations, see [49]. The remainder in (3.8) can be bounded, investigating the decay rate of the eigenvalues [26,49]. We obtain with C > 0, see [26, Proposition 2.5] and hence, An analytic covariance yields an exponential decay.
Except for several simple covariance functions, e.g., the exponential kernel, the eigenvalue problem (3.6) has to be solved numerically. This has been addressed, e.g., in [49,42]. To assure global differentiability, we employ a Galerkin approximation based on a B-spline space. This leads to a generalized discrete eigenvalue problem. Let I min := min Ī , I max := max Ī , and q, N ∈ N and consider a (quasi-uniform) sequence referred to as mesh τ N . To τ N we associate the B-spline space S q,k N of polynomials of degree q on each sub-interval of τ N with global continuity k ∈ N. For the standard iterative construction procedure of B-splines, see, e.g., [24]. The Galerkin approximation of (3.6) reads, following [49, p.111 for all v N ∈ S q,k N . By means of the numerically computed eigenpairs, the discrete Karhunen-Loève expansion reads Based on the approximation properties of B-splines, see, e.g., [30], the Galerkin error contribution can be bounded, following [49], as (3.14) f , with C M > 0 and q + 1 < l.

Application to the Nonlinear Magnetic Material Law.
In the following we employ a Karhunen-Loève expansion to discretize the stochastic magnetic material law. We thereby focus on the truncation error, i.e., assume that N is sufficiently large. As in the deterministic case we setf (i) := ∂ i sf (·, s). To ensure the solvability of the PDE the L 2 (Ω × I)-convergence off tof M is not sufficient and we additionally require thatf Proof. Based on Assumption 3 it follows by [49, Theorem 2.24] and (3.9), that for λ n = 0, a.s., with 0 < δ < 1/2 − 1/l and a positive constant C δ additionally depending on the covariance and |Γ|. Hence, the sum on the right-hand-side converges for all M . Moreover, we can choose M ≥ M 0 large enough such that for almost all ω ∈ Ω, s ∈ I, Let f M be a continuously differentiable prolongation off M from I to R + 0 , such that Assumption 2 is satisfied. This can be achieved, e.g., by the construction given in [31]. Then ν M (ω, s) := f M (ω, s) /s can be used as material coefficient for the stochastic problem. Complementary to this result, formulas assuring that the shape constraint (3.17) is verified, will be derived in the next section for a concrete setting.
Remark 2. We observe that we cannot model f or ν, as normal (or log-normal) random field, as their derivatives would not be bounded uniformly. In particular the trajectories of f would be non-monotonic with a probability greater than zero.

Practical Realization and Numerical Example.
In this section a practical realization under minimal assumptions on f is discussed. Measured data is supposed to be available at equidistant points I min ≤ŝ 1 <ŝ 2 < · · · <ŝ R ≤ I max . We introduce the table and assume that the data is monotonic, i.e.,f i1j ≤f i2j , for i 1 ≤ i 2 and j = 1, . . . , Q. The data is interpolated using C 1 monotonicity-preserving cubic splines, see [27]. For data with increased oscillations due to measurements a procedure as outlined in [44,47] should be used. The situation Q = 1 is very common in practice, and additional assumptions on the covariance are needed in this case. To solve the Karhunen-Loève eigenvalue problem, as outlined in Section 3.1, the correlation function Covf (s, s) Covf (t, t) Expected value E f and correlation function k f for the data given in [45].
is chosen to be approximated by the Gaussian kernel where L denotes the correlation length and δ > 0 a parameter. By rescaling with the we obtain the covariance. The associated discrete eigenpairs (λ n , b n ) N n=1 are obtained by discretization with S 3,1 N . Then, the random material relation is approximated by In (3.21) Ef is obtained by projecting the interpolated sample mean Ef N . Here, we determine M such that the relative information content satisfies For a more rigorous approach to the recovery of the Karhunen-Loève approximation from measured data by means of an a posteriori error analysis, we refer to [5]. Sample realizations of the Y n can be determined by (3.7), however here, we model them to be distributed uniformly as U − √ 3, √ 3 . The parameter δ is used to assure that (3.21) satisfies the shape constraints (3.17), see also [8, p. 1282] for a related discussion in the context of a linear material coefficient. In particular we only need to assure thatf M is monotonic. If Ef can be represented by a spline function, as is the case in the present setting, a simple condition for δ can be derived to this end. To simplify notation let η M (s) = denote the coefficient vector of Ef ∈ S 3,1 N . Then from [24] we recall that a sufficient condition for a B-spline to be monotonic is that its coefficients are increasing and hence, monotonicity can be assured by where we minimize only over those i with nonzero denominator. In a more general setting one could derive a similar expression by minimizing δ s given by over all s ∈Ī. Let us consider the material uncertainty of an electrical machine. In [45], measured data 1 , representing the material properties from twenty-eight machine stator samples (Q = 28) from the same production chain was presented. The interval of interest is given byĪ = [1, 1.55] and measurements were taken at R = 14 equidistant points. For the given criteria we truncate the Karhunen-Loève expansion with M = 3. Figure 2 depicts both the expected value and the correlation function. For illustration we compare the setting for three different correlation lengths L = 1/20, L = 1/10 and L = 1/2, respectively. Figure 3 depicts ten sample realizations for each correlation length, where the coefficient δ is chosen according to (3.23) and a uniform mesh with N = 60 spline basis functions is used. It can be readily observed, that a smaller correlation length, corresponding to trajectories with increased oscillations, demands for a smaller δ, i.e., a smaller perturbation magnitude. Also, the largest correlation length L = 1/2 gives the best agreement with the measured data and will be chosen in what follows.

Truncation Error and High-Dimensional Deterministic Problem.
We are now going to investigate the modelling error arising from a finite dimensional noise approximation, i.e., when ν is replaced by ν M . So far we have explained how this approximation can be achieved by means of the Karhunen-Loève expansion, however, hereafter, we do not restrict ourselves to this specific case anymore. For given ν M , for r ∈ R 3 , let h M (·, r) := ν M (·, |r|) r denote the associated vector function. Then A M is defined a.s. as the solution of Moreover, let f as well as f M satisfy Assumption 2, with constants α, α 0 , respectively. Then we have a.s.
Proof. As f satisfies Assumption 2, we have a uniform strong monotonicity property, i.e., a.s.
and because of equations (3.1) and (3.25) and the Cauchy-Schwarz inequality For the right-hand-side we further obtain Due to (3.26) we have control of the truncation error. For simplicity, this error is omitted in the following, i.e., we assume that the uncertain input has a finite dimensional noise representation: Assumption 4. The random field ν, resp. f , depends (continuously) on M independent random variables solely, i.e., a.s.
We recall that Y may also refer to random variables in closed-form representations of ν or coefficients in spline models, among others. Based on Proposition 3.2 and Assumption 4 it follows from the Doob-Dynkin Lemma [46] (cf. [4]) that we can write We recall that the Y n have a bounded image and a joint probability density function We now introduce the following assumption for f : Γ × R + 0 → R + 0 . Assumption 5. The stochastic material law f (y, ·) satisfies Assumption 1 for ρ-almost all y ∈ Γ, with constants α and β independent of y. The stochastic problem can be recast into a deterministic one with 3 + M dimensions.

A Stochastic Collocation
Method for the Nonlinear curl-curl Formulation. In the following, the variables x ∈ D and y ∈ Γ will be referred to as deterministic and stochastic variable, respectively. The solution of (3.37) requires discretization in both variables as well as a linearization procedure. To this end, we carry out: (i) deterministic discretization based on lowest order H (curl)-conforming finite elements with maximum stepsize h, (ii) stochastic discretization based on a collocation procedure on a tensor grid or sparse grid of level q, (iii) l-times iteration of the linearized system of equations by means of the Kacanov or Newton-Raphson method. As we will see in Section 4.4, the use of global, higher order polynomials over Γ is justified by the regularity of the solution, whereas the collocation procedure is particularly attractive for nonlinear problems due to the ease of implementation. We will then proceed by analyzing the approximation error originating from finite element discretization ε h , stochastic collocation ε q and linearization ε l , respectively. By the triangle inequality these errors can be decomposed as Additional sources of error can be identified, in particular quadrature errors and the error from numerically solving linear systems of equations. However, these errors will be omitted here. We also claim that all three steps of the proposed scheme commute. This has been shown for the deterministic case in [51] and generalizing to the stochastic collocation method is straightforward.

Galerkin Finite Element
Approximation. Equation (3.37) is approximated in the deterministic variable by the Galerkin finite element method. Higher order schemes are well established and could be employed, however, as our focus lies on the stochastic part, we restrict ourselves to lowest order schemes. We consider discretizations of the Lipschitz polyhedron D with a simplicial mesh T h , with maximum size h > 0. The mesh is assumed to be quasi-uniform in the sense of [14,Definition 4.4.13], in particular where B T denotes the largest ball contained in T . We then introduce the discrete spaces i.e., V h and W h are spanned by lowest order Nédélec and Lagrange elements, respectively. As in the continuous case, the space of (discrete) divergence free functions, is introduced as We observe thatV h is not a subspace ofV , as (weak) discrete divergence free functions are not (weak) divergence free in general. The deterministic finite element approximation consists in computing A h : Γ →V h such that for ρ-almost all y ∈ Γ, holds. Existence and uniqueness can be established based on a discrete Poincaré-Friedrichs inequality [32,Theorem 4.7], see, e.g., [53].

Stochastic Collocation Method.
Stochastic discretization is based on a collocation approach using either a tensor or a sparse grid, see, e.g., [52,6,7,38]. Starting with the (isotropic) tensor grid, following [38], the collocation points are given as where in each dimension m = 1, . . . , M , we have n(q) = p(q)+1 collocation points and N q = n(q) M in total. Note that p refers to the underlying polynomial degree, which we identify with the level for the tensor grid case as p(q) = q. Also, a global index k is associated to the local indices in the usual way [7]. The collocation points are chosen as the roots of the orthogonal polynomials associated to the probability density function ρ. As commonly done [39,7] we introduce the notation y = (y m ,ŷ m ) ,ŷ m = (y 1 , . . . , y m−1 , y m+1 , . . . , y M ). Let Q p (Γ m ) be the space of polynomials of degree at most p in Γ m . Then we introduce in each dimension the one-dimensional Lagrange interpolation operator I m p : where l i m (y m ) is the Lagrange polynomial of degree p associated to the point y i m . The tensor grid interpolation formula reads as (4.9) where l k (y) is the global Lagrange polynomial associated to the point y k ∈ H T q,M . An isotropic tensor grid, with n points in each direction, can only be used for moderate dimensions M , as the total number of collocation points grows as n M . Therefore, collocation in higher dimensions is based on sparse grids [13,41]. For simplicity we consider isotropic Smolyak grids, solely. Anisotropic sparse grids are discussed, e.g., in [40]. Following [38], let j ∈ N M 0 be a multi-index and now also depend on the multi-index j. For the choice p(j) = 2 j for j > 0 and p(0) = 0, the Smolyak formula is given by The associated sparse grid is denoted H S q,M . Evaluating (4.9) and (4.11) requires solving for all collocation points y k in H T q,M and H S q,M , respectively. 4.3. Linearization. At each collocation point, iterative linearization is carried out until the linearization error is found to be sufficiently small. The l-th iterate, l ∈ N, consists in computing A h,q,l ∈ Q q (Γ) ⊗V h such that for k = 1, . . . , N q where Q q (Γ) refers to the polynomial space associated either to tensor or to Smolyak interpolation. For a precise definition of these spaces, see, e.g., [7]. The representation (4.13) follows [25,16] and in particular we consider the linearization, (4.14) h L (·, r) = ν (·, |r l−1 |) r, for r, r l−1 ∈ R 3 . This is usually referred to as Kacanov method or successive substitution in the literature. Note that the case of the Newton-Raphson method, i.e., (4.15) h L (·, r) = ν (·, |r l−1 |) r + ν (1) (·, |r l−1 |) |r l−1 | r l−1 ⊗ r l−1 (r − r l−1 ) , is also covered. Under restrictions on the starting point A h,q,0 and damping, if necessary, A h,q,l converges to A h,q . At each step equation (4.13) is well-posed by the Lax-Milgram Lemma for both choices. For the Kacanov method, this follows by observing that For the Newton-Raphson method we refer to Lemma 4.1 below.

Stochastic Regularity and Convergence Analysis.
Convergence of the stochastic collocation method introduced above, can be established once the regularity of the solution is known. Whereas for the present nonlinear elliptic problem, the regularity w.r.t. the deterministic variable x is well known, to our knowledge, the stochastic regularity of the solution A w.r.t. y has not been investigated. In the case of a linear elliptic PDE it is well known, that under some mild assumptions the solution is an analytic function of the stochastic variable [7,39,23]. Similar results hold true for several types of nonlinear problems, see [19,20]. Here, the mapping (4.16) y → ν (·, | curl A (y) |) is real, but not complex differentiable, see [33,10], and this impedes a complex analysis. Moreover, the techniques presented in [19], based on the implicit function theorem, cannot be applied as a norm gap arises. More precisely, the nonlinearity h : L p (D) 3 → L q (D) 3 can be differentiated only for q < p, see [50]. Higher order differentiability even requires a larger difference between p and q. Therefore, we conduct an explicit higher order sensitivity analysis to precisely determine the stochastic regularity.
Proof. This result has been established, e.g., in [34, Lemma 3.1]. In particular this implies, that the bilinear form b d (u; ·, ·), defined by is continuous and coercive onV . We recall thatV h is not a subspace ofV . However, b d is continuous and coercive onV h , too. The form b d arises naturally when sensitivities are computed and also in the linearized system obtained by the Newton-Raphson method, which is well-posed at each iteration step by the properties just established. Provided that f is k-times differentiable, the question arises whether higher order derivatives of the solution exist as well. In a first step, we consider the case k ≤ 3. We impose the following assumption on the material law.
Assumption 6. For the parametric material law there holds with bounded and uniformly continuous derivatives, and additionally The vanishing higher order derivatives of f around the origin are used here to ensure differentiability in presence of the absolute value. Under the previous assumption we infer that D k r h is bounded. Lemma 4.2. Let Assumption 6 hold true. Then for |α| 1 ≤ k ≤ 3, ∂ α r h j is continuous and |∂ α r h j | ≤ C k . Proof. For r = 0, as h j (·, r) = f (·, |r|)r j /|r| we see that h j is k-times continuously differentiable. Hence, if ∂ α r h j (·, r) is bounded for r → ∞ and r → 0 the result follows. For r → ∞ we observe that f (k) is bounded (by Assumption 6) and that the same holds true for ∂ α r (r/|r|). For r → 0 we first observe that by the rule of l'Hôpital (4.21) ν (k−1) (·, 0) = f (k) (·, 0)/k and hence ν (1) (·, 0) = ν (2) (·, 0) = 0 by Assumption 6. Hence, using the expressions for D k r h given in Appendix A we infer that D 1 r h(·, 0)(s 1 ) = ν(·, 0)s 1 and that D 2 r h(·, 0) = D 3 r h(·, 0) = 0. By examining the proof of Lemma 4.2 we observe that if In the continuous case, formally differentiating the strong form of the boundary value problem, we obtain for the derivative A γ where, as shown in the Appendix, F k is given by r h α (·, curl A) curl A π1 , . . . , curl A π ca(π) .
With Π * we denote the set of partitions of γ − α, such that ca (π) > 1 if α = 0, where ca (π) refers to the cardinality of π. We refer to Appendix B for a more detailed definition of the underlying sets of F k . We observe, that π i < γ for i = 1, . . . , ca (π) and hence the derivatives contained in the right-hand-side are of lower order. Equation (4.22) is the basis for establishing the regularity of the solution with respect to the stochastic variable. This in turn determines the convergence rate of the corresponding discretization error. Before bounding the collocation error, we discuss the finite element error. To this end, we assume the following: Assumption 7. The solution of (3.37) has the additional regularity A ∈ L ∞ (Γ, H s (curl; D)), with s ∈ (1/2, 1], and the same holds true for the weak solution of (4.22), if it exists. Based on this regularity assumption, we can establish the following result: Lemma 4.3. Let Assumptions 5 and 7 hold true, then the deterministic error is bounded as where C depends on A and on s but is independent of y.
Proof. For the deterministic error ε h , also in the present nonlinear case, Céa's Lemma holds for ρ-almost all y ∈ Γ, see [10]. Then from [37, Theorem 5.41] we obtain where Π h is the canonical interpolation operator [37, p. 134] and the constant C 2 depends on s, but is independent of y. As A ∈ L ∞ (Γ, H s (curl; D)), we conclude that ε h ≤ C 3 h s . Based on the sensitivity analysis carried out in this section, we can now use Theorems 4 and 5 of [38] to establish the main result.
Theorem 4.4. Let Assumptions 5, 6 and 7 hold true. There holds where C is the constant from Lemma 4.3. Moreover, let i s be an integer 1 ≤ i s ≤ k, such that i s = 1 for s ∈ (1/2, 3/4), i s = 2 for s ∈ [3/4, 1) and i s = 3 for s = 1, respectively. For the isotropic tensor grid collocation method we have with respect to the level (polynomial degree) q and the number of collocation points N q , respectively. For the sparse grid collocation method and M ≤ i s , there holds with respect to the level q and the number of collocation points N q in the sparse grid, respectively. The constants C T , C S depend on s, A, ν, ρ and ξ ≈ 2.1.
Proof. The deterministic error estimate has been established in Lemma 4.3. In a first step, we bound the stochastic collocation error for the isotropic tensor grid. In this case, the collocation error can be recast as an interpolation error We will consider the case M = 1, solely, as the result for M > 1 follows by induction, see, e.g., [15,Lemma 7.1] or [38,Theorem 4]. The collocation error is related to the best-approximation error in Q q (Γ) ⊗V h , following [7], as The error decay then depends on the smoothness of A (and hence A h ) with respect to the stochastic variable. We note that, as M = 1, ∂ γ y A simplifies to ∂ is y A. Using (4.22) and the coercivity of b d we can formally bound
This follows from the Sobolev embedding theorem, see, e.g. [37, Theorem 3.7], as we have ∂ j y curl A(y) ∈ H s (D) 3 . (iii) For s = 1 we obtain i s = 3, as we can show that ∂ j y curl A(y) ∈ L 6 (D) 3 , for j = 0, 1, 2 using again the Sobolev embedding theorem. Hence, the results of Jackson quoted from [38] yield with C 4 depending on s and by induction for M > 1 The results for the sparse grid collocation error can be inferred from [38] or [13] once bounds on mixed derivatives of order k, i.e., ∂ γ y A with γ j ≤ k for j = 1, . . . , M have been established. Using the same arguments as above we obtain for |γ| 1 ≤ i s . This in turn ensures bounded mixed derivatives of order k = i s /M for the case M ≤ i s , solely. Then the results follow from Theorem 5 of [38]. Remark 3. Concerning the stochastic discretization error, in the two-dimensional case we can choose i s = 2 for s ∈ [1/2, 2/3) and i s = 3 for s ∈ [2/3, 1], respectively.
Remark 4. The linearization error could be included into this convergence estimate: provided that the initial values A h,q,0 (y k ) are sufficiently close to A h,q (y k ), there exists r ∈ (0, 1), such that (4.36) ε l ≤ Cr l for the linearization error of the solution A h,q,l of (4.13) obtained by the Kacanov method. An improved estimate could be obtained for the Newton-Raphson method.
Remark 5. The sparse grid convergence rate of O(N −γ is/M q ), with γ ∈ (0, 1), is smaller than the tensor grid convergence rate of O(N −is/M q ). This reduction, due to a lack of mixed regularity of the solution, is also observed in the numerical experiments in Section 5. However, for smooth solutions the sparse grid approach is expected to be more efficient for large M .
We now address the question whether a fast convergence, e.g., a rate of q −k for the tensor grid and arbitrary k ∈ N, can be obtained under suitable regularity assumptions on the material input data. This is true if curl A is bounded uniformly, e.g., for smooth domains and data, as can be seen by the preceding arguments. Also, if we accept a non-uniform constant with respect to the mesh size h, the decay of the stochastic discretization error can be improved, as stated in the following.
Proposition 4.5. Let Assumption 5 hold true and let ∂ α r ∂ β y h be continuous and bounded, for |α| 1 + |β| 1 ≤ k ∈ N. Then we have (4.37) ε q ≤ C T,h,k q −k , for a tensor grid and if M ≤ k, for a sparse grid, respectively. The constants additionally depends on A, ν, ρ.
Proof. As in the proof of Theorem 4.4 we bound F k as where we have used by the shape-uniformity of the mesh (4.2). This yields with |γ| 1 = k, for the solution A h (y) of (4.6). The result can now be established along the lines of the proof of Theorem 4.4 by observing that (4.33) holds true for k arbitrary.
A similar idea was applied in [43] in a two-dimensional setting to establish the Lipschitz continuity of the Newton operator.

Numerical Examples.
Several numerical examples are presented here to illustrate the findings. In 2D we consider a stochastic version of the p-Laplacian, where the solution is known exactly, and the L-shaped domain, respectively. In 3D the magnetostatic TEAM benchmark problem 13 [2] will be discussed with the stochastic material law of Section 3.2.1. Here, numerical results are obtained using FEniCS [51], and all non-uniform meshes are generated by Gmsh [28]. Tensor and sparse grids are generated using the Sparse grids Matlab kit [11,1].
In 2D, the curl operator reduces to curl (u e 3 ) = (∂ x2 u, −∂ x1 u, 0), where e 3 denotes the unit vector in the third dimension. The 2D equivalent of (2.9) then reads, find u ∈ H 1 0 (D) subject to Equation (5.1) is approximated by means of W h , i.e., lowest order nodal finite elements.

p-Laplace.
We consider D = (0, 1) × (0, 1), a constant current J = 2 and for s ∈ R + giving rise to the p-Laplace problem considered in [25]. The solution is given by as we use the inhomogeneous Dirichlet boundary condition u| ∂D . Modeling Y = p as a single random parameter, with Y > 1 a.s., we obtain a stochastic problem. For x 1 x 2  each Y we have u(Y, ·) ∈ H 1 (D) and even grad u(Y, ·) ∈ L ∞ (D) 2 . Hence, we expect a fast spectral stochastic convergence. Note that (5.2) violates the assumptions on the reluctivity for s → {0, ∞}. However, monotonicity and continuity results can be obtained as outlined in [17]. We model Y as uniformly distributed on (3,5), i.e., Y ∼ U (3, 5). Using a uniform triangulation, (5.1) is iteratively solved by means of the Kacanov method (with damping) until the solution increment is below 10 −12 in the discrete L ∞ -norm, which yields a negligible linearization error compared to the other sources of error. In Figure 4 we depict the expected value and variance of the solution, respectively. In Figure 5 the error E[u − u q,h ] H 1 0 (D) is depicted for different values of h. As we are not aware of a closed form solution of E[u], we approximate it by E q [u], with q = 10. From Figure 5 an exponential decay of the stochastic error can be observed until the corresponding discretization error level is attained. applied together with the constant excitation J = 10 5 . We adopt the material model

Var[u]
from [21]. Given initial values a 0 = 1.78, b 0 = 14, c 0 = 6000 and d 0 = 245, we set a = a 0 (1 + 0. Figure 6 on the left and right we depict the solution for Y 1 = Y 2 = 0 and random realizations of the reluctivity, respectively. Linearization is carried out as in the previous example. A coarse finite element (FE) grid as depicted in Figure 6 on the left, referred to as FE grid 1, is uniformly refined two (FE grid 2) and four times (FE grid 3), respectively. The stochastic error for a tensor grid H T q,2 , w.r.t. both the polynomial degree and the number of grid points, is shown in Figure  7. As a reference, a higher order polynomial approximation (q = 9) in the stochastic variable is used on each grid. For this example there holds grad u ∈ H s (D) 2 , with s = 2/3 − , with > 0. According to Theorem 4.4 in 2D, see Remark 3, at least a decay of q −2 is expected as s ∈ [1/2, 2/3). For the finest FE grid we numerically observe a decay of even q −2.987 . Hence, the convergence rate seems to be insensitive to the small parameter and the result could possibly be slightly improved. For grids 1 and 2 the decay is also faster than predicted. In Figure 8 the errors for a sparse H S q,2 are depicted. As i s /M = 1 in this case, the convergence is assured. We observe an algebraic decay w.r.t. to the number of grid points. Also the tensor grid approach is more efficient in this case as, e.g., for N q = 49 and FE grid 3, the error using H S q,2 is 4.12 × 10 −4 , whereas the error using H T q,2 is 1.77 × 10 −5 . 5.3. TEAM Benchmark. TEAM benchmarks are setup to validate electromagnetic codes and in particular magnetic field simulations. Here, we investigate the nonlinear magnetostatic TEAM 13 problem. A magnetic field in three thin iron sheets is generated by a rectangular coil, with blended corners, as depicted in Figure 9 on the right. Due to symmetry, only the upper half of the iron sheets is visualized. In a pre-processing step, an electrokinetic problem is solved to obtain the source current distribution J with a total imposed current of 3000A per cross section. Gauging is enforced through a Lagrange multiplier and a mixed formulation see, e.g., [36]. The computational domain is truncated, applying homogeneous Dirichlet boundary conditions at a boundary, sufficiently far away from the problem setup. Strictly speaking   this iron-air interface problem would require minor modifications of the analysis presented in this paper, as mentioned in Remark 1. We replaced the material properties of the original benchmark and employed the stochastic B-H curve presented in Section 3.2.1, with L = 1/2 and δ = 2 instead. Extrapolation beyond the data range is carried out as described in [47]. Note that for this case we are concerned with C 1 trajectories, solely. A tetrahedral mesh, see also Figure 9, is generated using the software Gmsh [28] and two steps of uniform refinement are carried out. We refer to the different FE grids with 7388, 41958 and 261196 total degrees of freedom as grid 1, 2 and 3, respectively. In Figure 9 on the left the magnetic flux density distribution is depicted for the expected value of the B-H curve. The nonlinear problem is linearized and iterated as explained in the p-Laplace example with a linearization increment below 10 −5 in the discrete L ∞ -norm. In view of the expected low regularity of the solution and the small number of random inputs a tensor grid H T q,3 is used. The stochastic convergence is depicted in Figure 10 for all three FE grids. Here, the stochastic error is estimated as H(curl)-norm of E[A h,q,l ] − E[A h,q+1,l ], as we do not have an analytical solution. The initially rapid convergence deteriorates until the linearization error level is attained. We also observe that the convergence is algebraic. However, the convergence is faster than the predicted rate q −1 in view of the limited differentiability of the reluctivity.   6. Conclusions. In this work, we have addressed the stochastic nonlinear elliptic curl-curl equation with uncertainties in the material law. Assumptions on the input have been formulated in order to obtain a well-posed stochastic formulation and it was shown that they can be fulfilled when a suitable discretization of the random input is carried out by the truncated Karhunen-Loève expansion. As monotonicity is required for the trajectories of the material law, oscillations can only occur with a rather small magnitude. A stability result for approximations of the random input was also derived. In the second part of the paper, a stochastic collocation method for the nonlinear curlcurl equation was analyzed. Under moderate differentiability assumptions on the material law, a convergence rate of q −k was obtained for the stochastic collocation error using tensor grids, where 1 ≤ k ≤ 3. For smooth boundaries and data this estimate holds true for all k ∈ N. These estimates were shown to be in good agreement with numerical results for academic and benchmark examples. Convergence results for sparse grids were also obtained. However, in this case convergence can be expected only for a limited number of random input parameters.