Laurent Baratchart

Sylvain Chevillard

Douglas Hardin

Juliette Leblond

Eduardo Andrade Lima

Jean-Paul Marmorat

Magnetic moment estimation and bounded extremal problems

Keywords: Inverse moment estimation problems, elliptic partial differential equations, best constrained approximation in Hilbert spaces, harmonic functions, geosciences and planetary sciences, paleomagnetism 35J15, 35R30, 35A35, 42B37, 86A22; Secondary: 45K05, 46F12, 47N20

We consider the inverse problem in magnetostatics for recovering the moment of a planar magnetization from measurements of the normal component of the magnetic field at a distance from the support. Such issues arise in studies of magnetic material in general and in paleomagnetism in particular. Assuming the magnetization is a measure with L 2 -density, we construct linear forms to be applied on the data in order to estimate the moment. These forms are obtained as solutions to certain extremal problems in Sobolev classes of functions, and their computation reduces to solving an elliptic differential-integral equation, for which synthetic numerical experiments are presented.

Introduction

In this paper, we study a family of constrained best approximation problems in certain Hilbert spaces of harmonic gradients, which are motivated by inverse magnetization problems arising in geosciences and planetary sciences. The physical setting, along with instrumentation issues and related questions from paleomagnetism, are described in the introductory sections of [START_REF] Baratchart | Silent and equivalent magnetic distributions on thin plates[END_REF][START_REF] Baratchart | Asymptotic method for estimating magnetic moments from field measurements on a planar grid[END_REF][START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF][START_REF] Lima | Fast inversion of magnetic field maps of unidirectional planar geological magnetization[END_REF] to which this work may be viewed as a sequel. The goal is to analyze magnetized thin rock samples carrying some unknown magnetization distribution to be estimated from measurements of the magnetic field. Specifically, a sample S is identified with a compact subset of a plane P 0 ⊂ R 3 , and is assumed to support a R 3 -valued magnetization distribution m. The measurements consist of pointwise values, on a compact set Q contained in a plane P h parallel to P 0 lying at distance h > 0 from it, of the normal component of the magnetic field produced by m (normal to P h , that is). For definiteness, we assume that both S and Q are compact rectangles, centered at the origin, in the horizontal planes R 2 × {0} ⊂ R 3 and R 2 × {h}, respectively. Thus, the measured quantity is the vertical component b 3 [m] of the magnetic field produced by m, evaluated on Q. Maxwell's equations for the magnetostatic case [START_REF] Jackson | Classical Electrodynamics[END_REF]Ch. 5] imply that this magnetic field derives from a scalar magnetic potential Λ which is a solution to the Poisson equation ΔΛ = div m in R 3 . Consequently, we have that b 3 [m] = -µ 0 ∂ x 3 Λ with µ 0 = 4π × 10 -7 , see e.g. [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF][START_REF] Baratchart | Silent and equivalent magnetic distributions on thin plates[END_REF].

It is well-known that some magnetizations generate no external magnetic field (i.e., the zero field): these are called silent. A description of silent magnetization distributions supported in a plane is given in [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF], under weak assumptions on the distribution. In particular, those having compact support admit a simple characterization, cf. Section 2.2. Moreover, it is shown in [START_REF] Baratchart | Silent and equivalent magnetic distributions on thin plates[END_REF] that if the normal component of the field vanishes everywhere on an open subset of a plane, then the field is identically zero and therefore the magnetization is silent. Thus, with the previous notation, it holds that b 3 [m] = 0 on Q if and only if m is a silent magnetization of the type described in [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF]. The existence of such magnetizations is a major obstacle for solving the inverse magnetization problem, for it implies that additional assumptions must be made to achieve full recovery. In the present work, we focus on a much simpler inverse problem, which consists in estimating the net moment m of the unknown magnetization m, i.e. the integral1 of m. So, we merely aim at recovering a vector in R 3 (namely m ) rather than the whole distribution m on the rectangle S. Unlike magnetization recovery, net moment recovery is meaningful without extra-assumptions on m, because silent magnetizations are easily seen to have zero moment.

Though much less ambitious than solving the full inverse magnetization problem, determining the net moment is still a significant goal. Indeed, classical magnetometers estimate the net moment of a sample by comparing it to a dipolar source, an approximation which is valid at a distance from the sample only. When dealing with weakly magnetized objects, the field away from the sample gets too small and is easily dominated by other, spurious magnetic sources, so that measurements by standard magnetometers may become unreliable. This is a typical issue in paleomagnetism, which justifies our present attempt to devise algorithms to estimate the moment from measurements of (the vertical component of) the field, close to the sample. Such measurements can be performed by very sensitive instruments known as Scanning QUantum Interference Device (SQUID) microscopes, see [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF][START_REF] Lima | Fast inversion of magnetic field maps of unidirectional planar geological magnetization[END_REF] for more explanation.

From a physical viewpoint, knowing the net moment of a magnetization m is a basic piece of information that is relevant in itself. It gains further significance when putting additional assumptions on m, assuming for instance that it is unidimensional (i.e. that it takes values in a 1-dimensional subspace of R 3 ). Unidimensionality is a fairly common assumption in paleomagnetism, at least for igneous rocks that cooled down in the presence of an ambient field during formation and whose remanent magnetization keeps record of the strength and direction of that field. Clearly, if a magnetization is unidimensional then its direction must be the one of the net moment. Now, if we are able to compute this moment, one may consider rather efficient regularizing techniques in the Fourier domain to recover the magnetization, see [START_REF] Lima | Fast inversion of magnetic field maps of unidirectional planar geological magnetization[END_REF] and [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF]Sec. 4.1]. More can be said in this connection: in fact, to any magnetization m compactly supported in a plane P 0 , there is a unidimensional magnetization u supported on P 0 which generates the same field on a given connected component of the complement of P 0 . Surprisingly perhaps, the direction of u can even be picked arbitrarily provided it is not tangent to the plane [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF]Thm 3.6]; however, clearly, u will not have the same support as m in general, and will typically be non-compactly supported. Now, even in this case, choosing u to have the direction of m has a regularizing effect on the Fourier techniques just mentioned, because the singularity of the Fourier transform û at 0 can then be singled out explicitly, compare [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF]Eqns ( 29)-( 30))]. This dangles the prospect of efficiently computing a unidimensional magnetization which differs from the true one by a (generally non-compactly supported) silent magnetization. Because the latter can be characterized completely, estimating the net moment may thus be envisaged as an initial step to study the inverse magnetization problem in general.

Having made the case for net moment recovery and pointed out that it was an easier task than full inversion, we should nevertheless stress that it is nontrivial. A main reason is that measurements of b 3 [m] are taken merely on Q which does not surround S. If measurements were available on the entire plane R 2 × {h} (which does surround S if one takes into account points at infinity), then asymptotic formulas from [START_REF] Baratchart | Asymptotic method for estimating magnetic moments from field measurements on a planar grid[END_REF] would allow us in principle to approximate the net moment arbitrary well, at least if m is smooth enough, say if it is a measure. However, the fact that measurements are available on Q only makes the problem ill-posed, in that small differences in b 3 [m] on Q may result in large differences on m , cf. Section 2.2. This is precisely the issue that we address in the present work. Further reasons why net moment recovery is difficult in practice lie with the presence of noise in the measurements, inherent in such problems. In this paper, we only deal with synthetic examples corrupted by a small additive Gaussian white noise, and we do not touch upon this important facet of the situation except for some general comments in Sections 2.3 and 5.2.

We shall restrict ourselves to the elementary case where m is a measure supported on S with square summable density there. This makes for a Hilbertian framework which keeps the analysis simple. Hopefully this case is already typical of the main features of our method, though there is grounds to develop a similar approach in non Hilbertian contexts, for instance when magnetizations are measures normed with the total variation. Roughly speaking, what we do in this paper is to construct a R 3 -valued function

φ on Q such that Q b 3 [m] φ ∼ S m = m (1)
for all m ∈ L 2 (S, R 3 ) of given norm, with a bound on the norm of the derivative of φ in L 2 (Q, R 3×2 ) which serves as a regularization parameter. In other words, we construct a linear estimator φ to be applied to the data b 3 [m] for estimating m , but a trade-off exists between quality of approximation in (1) and oscillation of φ. Note that oscillatory behavior is undesirable to evaluate the left hand side of (1) from pointwise values of b 3 [m] with good accuracy. Besides, our choice to bound the derivative of φ makes for an easy theoretical solution of the extremal problems involved in our analysis.

To conclude this introduction, let us mention also that when m is the sum of several magnetization distributions with disjoint supports, then the individual moment of a component can be estimated even though we measure the global field, provided we have a priori knowledge of a neighborhood of the support of that component that does not intersect the support of any other. This can be done by essentially the same techniques as those we use to get [START_REF] Adams | Sobolev Spaces[END_REF]. Computing such local moments could be of interest to get information on the support of m, for instance to confirm that a region where no magnetization is expected has indeed zero or near-zero moment.

The outline of this work is as follows. Section 2 is devoted to notation, definitions and preliminaries, including a description of the net magnetic moment recovery problem. This issue is recast in Section 3 as a bounded extremal problem for which existence and uniqueness results are proven. A characterization of the solution in terms of a critical point equation involving a Lagrange parameter is then derived. Two methods for solving the critical point equation numerically are presented in Section 4, one of which is used in Section 5 to provide numerical illustrations. Section 6 has concluding remarks.

Notation and preliminaries

Function spaces

For x = (x 1 , • • • , x n ) t ∈ R n , we let x = (x 2 1 + • • • + x 2 n )
1 2 be its Euclidean norm; here and below, a superscript "t" means "transpose". Denote the scalar product of x and y by x.y = j x j y j . The notation stands irrespective of n but no confusion should arise, and we are only concerned with n = 2, 3. We write E for the closure of a set E and for n = 2 we indicate the open disk of center x and radius R with D(x, R).

We write C ∞ (Ω) for the space of smooth functions on an open set Ω ⊂ R n having continuous derivatives of any order, and C ∞ 0 (Ω) for the space of smooth functions with compact support in Ω. Recall that distributions on Ω are linear forms on C ∞ 0 (Ω) which are continuous for a certain topology, the precise definition of which may be found in [START_REF] Schwartz | Théorie des distributions[END_REF]Sec. I.2]. The support of a distribution d on Ω, denoted by supp d, is the largest closed set E ⊂ Ω such that d(ϕ) = 0 for all ϕ ∈ C ∞ 0 (Ω \ E). Distributions get differentiated by the usual rule:

∂ x i d(ϕ) = -d(∂ x i ϕ), with ∂ x i to
mean the partial derivative with respect to the coordinate x i . It is standard to regard a locally integrable function ψ as a distribution by putting ψ(ϕ) = ψϕ. For φ a distribution or a function, we let ∇φ = (∂ x 1 φ, • • • , ∂ xn φ) t for the gradient of φ. The divergence operator, which is the distributional adjoint to ∇, is denoted as ∇. and acts on R n -valued distributions as

∇.(φ 1 , • • • , φ n ) = Σ 1≤j≤n ∂ x j φ j .The Laplace operator Δ is defined by Δ φ = Σ 1≤j≤n ∂ 2
x j φ. These operator notations stand irrespective of n, but no confusion should arise. In fact we use them for n = 2, except in Equation (4) when we introduce the magnetic potential in general form, and after Equation (10) when we introduce silent magnetizations.

We need some standard functions spaces but only in dimension two, so we restrict the discussion to this case. For Ω ⊆ R 2 , we denote by L 2 (Ω, R k ) the familiar Lebesgue space of R k -valued square summable functions which is a Hilbert space with norm and scalar product:

φ 2 L 2 (Ω,R k ) := Ω φ 2 d , φ, ψ L 2 (Ω,R k ) := Ω φ.ψ d , with d to indicate Lebesgue measure. We simply write L 2 (Ω) if k = 1.
If Ω is open, we put W 1,2 (Ω) for the standard Sobolev space of functions belonging to L 2 (Ω) together with their first distributional derivatives. It is a Hilbert space with norm:

φ 2 W 1,2 (Ω) = φ 2 L 2 (Ω) + ∇φ 2 L 2 (Ω,R 2 ) . We let W 1,2 0 (Ω) be the closure of C ∞ 0 (Ω) in W 1,2 (Ω).
If Ω is bounded, the Sobolev-Poincaré inequality [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]Cor. 9.19] implies there exists a constant C > 0 (depending on Ω) such that

φ L 2 (Ω) ≤ C ∇ φ L 2 (Ω,R 2 ) , ∀φ ∈ W 1,2 0 (Ω) . ( 2
)
Therefore

• W 1,2 (Ω) and ∇ • L 2 (Ω,R 2 )
are equivalent norms on W 1,2 0 (Ω), namely:

∇ u L 2 (Ω,R 2 ) ≤ u W 1,2 (Ω) ≤ √ C 2 + 1 ∇ u L 2 (Ω,R 2 ) , u ∈ W 1,2 0 (Ω). ( 3
)
When Ω is bounded and Lipschitz-smooth, meaning that its boundary ∂Ω is locally the graph of a Lipschitz function, then functions in W 1,2 (Ω) have a trace on ∂Ω which belongs to the fractional Sobolev space W 1/2,2 (∂Ω), consisting of square summable functions g with respect to arclength on ∂Ω for which

g W 1/2,2 (∂Ω) = g L 2 (∂Ω) + ∂Ω×∂Ω |g(t) -g(t )| 2 t -t 2 |dt| |dt |) 1/2 < +∞,
where |dt| indicates (the differential of) arclength. The trace operator is continuous and surjective onto W 1/2,2 (∂Ω). Functions with zero trace are exactly those belonging to W 1,2 0 (Ω), and they are also those whose extension by zero outside Ω defines a function in W 1,2 (R 2 ). We also make use of the fractional Sobolev space W 3/2,2 (Ω). The latter consists of functions g ∈ W 1,2 (Ω), each partial derivative of which satisfies:

Ω×Ω |∂ x i g(x) -∂ x i g(y)| 2 x -y 3 dxdy) 1/2 < +∞, i = 1, 2.
The Sobolev space W 2,2 (Ω), consists of functions in L 2 (Ω) whose partial derivatives of the first order lie in W 1,2 (Ω). It is a Hilbert space with norm

φ 2 W 2,2 (Ω) = φ 2 L 2 (Ω) + ∇φ 2 L 2 (Ω,R 2 ) + 1≤i≤j≤2 ∂ x i ∂ x j φ 2 L 2 (Ω)
.

We refer to [START_REF] Adams | Sobolev Spaces[END_REF], [8, Ch. 9], [23, Ch. V, VI] for classical properties of Sobolev spaces. For 0 ≤ α < 1, we put C α (Ω) for the space of Hölder continuous functions of exponent α on Ω.

When f is a function on Ω, we denote by f the extension of f by zero outside Ω. That is, f is equal to f on Ω and to 0 elsewhere. We also use at places the symbol ∨ to mean concatenation with another function defined on R 2 \ Ω:

for instance f = f ∨ 0. Clearly, f L 2 (Ω) = f L 2 (R 2 ) and if f ∈ W 1,2 0 (Ω) then f ∈ W 1,2 (R 2 ).

Operators involved

The magnetic potential associated with a magnetization, modeled as a compactly supported R 3 -valued distribution m on R 3 , is given by

Λ(x) = - 1 4π ∇.m , 1 . -x = 1 4π m , ∇ 1 . -x , x ∈ R 3 \ supp m, (4)
where the brackets denote the duality pairing between distributions and C ∞ 0 (R 3 )-functions. Note that the formula makes sense because y → 1/ yx coincides with a C ∞ 0 (R 3 )function on a neighborhood of supp m when x / ∈ supp m. In the particular case where m is a R 3 -valued measure supported on a compact subset S of the horizontal plane R 2 × {0}, with density m = (m 1 , m 2 , m 3 ) there, then we get that

Λ(x) = 1 4π S m(t) . (x -t) x -t 3 d (t), x ∈ R 3 \ S. (5)
Let P x 3 be the Poisson kernel of the upper half-space {x 3 > 0}, see e.g. [24, Ch. I]:

P x 3 : x → x 3 2π x 3 . ( 6
)
On this half-space, the previous expression for Λ implies that it is half the convolution over R 2 of P x 3 with the function R 1 m 1 +R 2 m 2 + m 3 , where R 1 , R 2 are the Riesz transforms on R 2 , see [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF]Sec. 2]). In other words, for x 3 > 0, 2Λ(x 3 ) is the harmonic extension of R 1 m 1 + R 2 m 2 + m 3 to the upper half space. Likewise, for x 3 < 0, 2Λ(x 3 ) is the harmonic extension of R 1 m 1 + R 2 m 2 -m 3 to the lower half space. Hereafter, we let S and Q be two Lipschitz-smooth bounded non-empty open sets in R 2 × {0} and R 2 × {h} respectively. For simplicity, we also assume that S is simply connected. One may think of S, Q as being rectangles whose closures S, Q are the sets S, Q introduced in Section 1. Observe that expression [START_REF] Baratchart | Silent and equivalent magnetic distributions on thin plates[END_REF] for Λ may as well be computed as an integral over S because the boundary ∂S = S \ S has Lebesgue measure zero. We often identify S and Q with subsets of R 2 , while the third component (0 in the case of S and h in the case of Q) is treated as a parameter.

The main operator under consideration here is m → b 3 [m], which maps the magne-

tization m = (m 1 , m 2 , m 3 ) ∈ L 2 (S, R 3 ) to the vertical component of the magnetic field on Q. Specifically, b 3 : L 2 (S, R 3 ) → L 2 (Q) is defined by b 3 [m] = -µ 0 ∂ x 3 Λ, x 3 = h. ( 7
)
Although the target space is here L 2 (Q), it is clear since h > 0 that the range of b 3 consists of restrictions to Q of smooth (even analytic) functions on {x 3 > 0}. Basic properties of the operator b 3 were given in [START_REF] Baratchart | Silent and equivalent magnetic distributions on thin plates[END_REF]Sec. 3], [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF]. Below, we recall some of them that will be used in the course of the paper.

Recall the notation m i = m i ∨ 0 ∈ L 2 (R 2 ). From the relation

∂ x j P x 3 (x) = ∂ x j x 3 2π x 3 = ∂ x 3 x j 2π x 3 for j = 1, 2 , we get for all m = (m 1 , m 2 , m 3 ) ∈ L 2 (S, R 3 ) that b 3 [m] = - µ 0 2 (∂ x 1 P h ) m 1 + (∂ x 2 P h ) m 2 + (∂ x 3 P x 3 ) |x 3 =h m 3 |Q . ( 8
)
It is easily checked that b 3 is a bounded operator. In fact, being an integral operator whose kernel is uniformly continuous on

S × Q, it is compact [17, Ch. III, Ex. 4.1]. Its adjoint b * 3 : L 2 (Q) → L 2 (S, R 3 ) acts on φ ∈ L 2 (Q) by the formula b * 3 [φ] = µ 0 2     (∂ x 1 P h ) φ (∂ x 2 P h ) φ -(∂ x 3 P x 3 ) |x 3 =h φ     | S , (9)
and is likewise bounded, even compact because so is b 3 [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]Thm 4.10]. A majorization of the norm is [START_REF] Baratchart | Silent and equivalent magnetic distributions on thin plates[END_REF]Sec. 3.3]:

b * 3 ≤ b with b = µ 0 2 4 √ 2 3 3/2 h . ( 10
)
Moreover, by analyticity properties of the Poisson extension, it is not difficult to check that b * 3 is injective whence b 3 has dense range in L 2 (Q) [START_REF] Baratchart | Silent and equivalent magnetic distributions on thin plates[END_REF]Sec. 3.3]. It follows from [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF]Thm 2.3] that ∇Λ ≡ 0 on R 3 \ S if and only if m 3 = 0 and ( m 1 , m 2 ) t is divergence free on R 2 , in the distributional sense. Such magnetizations are called silent because they generate the zero field outside of their support. Moreover, by [START_REF] Baratchart | Silent and equivalent magnetic distributions on thin plates[END_REF]Prop. 2], it holds that b 3 [m] = 0 if and only if m ∨ 0 is silent. This means that m 3 = 0 and that (m 1 , m 2 ) t belongs to the space D S ⊂ L 2 (S, R 2 ) consisting of vector fields whose extension by zero outside of S defines a divergence-free vector field on R 2 . In the proof of [5, Prop. 2], it is shown that D S can be parameterized by Sobolev functions:

D S = (-∂ x 2 ψ, ∂ x 1 ψ) t , ψ ∈ W 1,2 0 (S) ⊂ L 2 (S, R 2 ). (11)
Besides, it follows at once from [START_REF] Dahlberg | Weighted norm inequalities for the Lusin area integral and the nontangential maximal functions for functions harmonic in a Lipschitz domain[END_REF] and the density of C ∞ 0 (Ω) in W 1,2 0 (Ω) that the orthogonal space D ⊥ S to D S in L 2 (S, R 2 ) is comprised of those vector fields satisfying the distributional Schwarz rule. Since S is simply connected, these are the distributional gradients on S, by [22, Ch. II, Sec. 6, Thm VI]. Thus, appealing to [START_REF] Demengel | Espaces fonctionnels, Utilisation dans la résolution des équations aux dérivées partielles[END_REF]Thm 6.74] which implies that a distribution with L 2 derivative on a Lipschitz open set must be a Sobolev function 2 , we conclude that D ⊥ S is the set ∇W 1,2 (S) of gradients of Sobolev functions. Altogether, this leads us to a characterization of the kernel of b 3 and of the range

of b * 3 in L 2 (S, R 3 ) as follows (compare [5, Lem. 4]). If we set D S = Ker b 3 , then      D S = (-∂ x 2 ψ, ∂ x 1 ψ, 0) t , ψ ∈ W 1,2 0 (S) ⊂ L 2 (S, R 3 ) D ⊥ S = Ran b * 3 = b * 3 W 1,2 0 (Q) = ∇ W 1,2 (S) × L 2 (S) ⊂ L 2 (S, R 3 ) , ( 12
)
where D ⊥ S stands for the orthogonal space to D S in L 2 (S, R 3 ) and the second equality in the second line of [START_REF] Demengel | Espaces fonctionnels, Utilisation dans la résolution des équations aux dérivées partielles[END_REF]

holds because W 1,2 0 (Q) is dense in L 2 (Q). Also, rotating pointwise by π/2 the relation D ⊥ S = ∇W 1,2 (S), we get that ∇ W 1,2 0 (Q) ⊥ = (-∂ x 2 ψ, ∂ x 1 ψ) t , ψ ∈ W 1,2 (Q) ⊂ L 2 (Q, R 2 ) , 2
The proof given there for bounded C 1 -smooth domains carries over to the Lipschitz case.

hence vector fields in ∇ W

Moment recovery issues

By definition, the net moment m i of the i-th component m i of m ∈ L 2 (S, R 3 ) is given, for i = 1, 2, 3, by the formula

m i = S m i (t) dt = m , e i L 2 (S,R 3 )
with e 1 = (χ S , 0, 0) t , e 2 = (0, χ S , 0) t and e 3 = (0, 0, χ S ) t , where χ S denotes the characteristic function of S. Note that L 2 (S, R 3 ) ⊂ L 1 (S, R 3 ) because S is bounded, hence the above formula for m i makes sense. The net moment m of m is simply the vector

( m 1 , m 2 , m 3 ) t ∈ R 3 . All magnetizations m such that b 3 [m ] = b 3 [m] have the same net moment as m. Indeed, b 3 [m] = 0 ⇔ m ∈ D S ,
and since e 1 , e 2 are gradients of Sobolev functions on S (namely of x → x 1 , x → x 2 ), we see from [START_REF] Demengel | Espaces fonctionnels, Utilisation dans la résolution des équations aux dérivées partielles[END_REF]

that each element of D S is orthogonal to e i in L 2 (S, R 3 ) for 1 ≤ i ≤ 3. Hence, m is uniquely determined by b 3 [m], in other words there is a linear map µ : Ran b 3 → R 3 such that µ(b 3 [m]) = m for m ∈ L 2 (S, R 3
). This map, however, is not continuous. Otherwise indeed, it would have a continuous extension L 2 (Q) → R 3 by the Hahn-Banach theorem and thus, to each i ∈ {1, 2, 3}, there would exist

φ i ∈ L 2 (Q) such that the quantity b 3 [m] , φ i L 2 (Q) -m i = b 3 [m] , φ i L 2 (Q) -m , e i L 2 (S,R 3 ) = m , b * 3 [φ i ] -e i L 2 (S,R 3 ) (13)
vanishes for all m ∈ L 2 (S, R 3 ). But the last term in [START_REF] Ern | Theory and Practice of Finite Elements[END_REF] implies that to each ε > 0 there is

φ i,ε ∈ L 2 (Q) such that b * 3 [φ i,ε ] -e i L 2 (S,R 3 ) ≤ ε.
Then, computing as in [START_REF] Ern | Theory and Practice of Finite Elements[END_REF], we get that

b 3 [m] , φ i,ε L 2 (Q) -m i = m , b * 3 [φ i,ε ] -e i L 2 (S,R 3 ) ≤ ε m L 2 (S,R 3 ) . ( 14
)
Equation ( 14) shows we can estimate m with arbitrary relative precision, at least in principle, by taking the scalar product of the data b 3 [m] with some appropriate estimator φ i,ε . However we necessarily have that φ i,ε L 2 (Q) → +∞ when ε → 0, otherwise taking a weakly convergent subsequence would imply in the limit that e i ∈ Ran b * 3 , a contradiction. Thus, when making ε small in order to improve accuracy of the estimate b 3 [m], φ i,ε L 2 (Q) for m i , we tend to increase the norm of the estimator φ i,ε which in turn amplifies the effect of measurement errors when computing the estimate. This is a familiar situation with inverse problems, which calls for a regularization technique to find a trade off between accuracy of the estimate and precision in the computation; this trade off stands analog, in a deterministic context, to the classical compromise between bias and variance from stochastic identification.

The regularization method we will use is somewhat dual to the Tychonov-type. Yet, instead of bounding the norm of the unknown magnetization, m, we will control the W 1,2 0 (Q)-norm of the estimator. The reason for this choice is threefold. First, it makes sense to keep the norm of the derivative of the estimator at a low level to prevent oscillations, since the latter may spoil the evaluation of the net moment from pointwise values of b * 3 [m]. Second, designing the estimator so as to vanish on the boundary ∂Q puts smaller weight on the values of the field close to ∂Q, which have poor signal/noise ratio because of its rapid decay when the distance to the sample increases. Last but not least, constraining the W 1,2 0 (Q)-norm of the estimator makes for a relatively simple approximation problem to solve, with a critical point equation of elliptic type showing interesting regularity properties.

Hereafter we fix e ∈ Ran b * 3 ⊂ L 2 (S, R 3 ). In connection with net moment estimation, e can be any of the e i for i ∈ {1, 2, 3} and then we are in the generic case where e ∈ Ran b * 3 , but perspective will be gained if we discuss the more general case. First, observe in view of ( 12) that inf

φ∈W 1,2 0 (Q) b * 3 [φ] -e L 2 (S,R 3 ) = 0 . Second, whenever φ n ∈ W 1,2 0 (Q) is such that b * 3 [φ n ] -e L 2 (S,R 3 ) → 0 as n → ∞, then ei- ther e ∈ b * 3 W 1,2 0 (Q) or ∇φ n L 2 (Q,R 2 ) → ∞. For if ∇φ n L 2 (Q,R 2 )
is bounded, extracting a subsequence converging weakly to some φ ∈ W 1,2 0 (Q) and using that the restriction map b

* 3 : W 1,2 0 (Q) → L 2 (S, R 3
) is a fortiori continuous, by the Sobolev-Poincaré inequality (2), gives us in the limit that e = b * 3 (φ). To recap, we see from ( 12) that the quantity

b 3 [m] , φ L 2 (Q) -m , e i L 2 (S,R 3 ) ≤ b * 3 [φ] -e i L 2 (S,R 3 ) m L 2 (S,R 3 ) (15)
can be made arbitrarily small for appropriate

φ, but if e / ∈ b * 3 [W 1,2 0 (Q)] this is at the cost of letting ∇ φ L 2 (Q,R 2 ) grow unbounded.
In practice, one really measures b 3 [m] + n where n is some noise (deterministic or stochastic) and thus, we compute b

3 [m] + n , φ L 2 (Q) rather than b 3 [m] , φ L 2 (Q) . If the noise is treated in a deterministic fashion, a relevant estimate is not (15) but b 3 [m] + n, φ L 2 (Q) -m, e i L 2 (S,R 3 ) ≤ b * 3 [φ] -e i L 2 (S,R 3 ) m L 2 (S,R 3 ) + n L 2 (Q) φ L 2 (Q) .
So, given a priori bounds on m L 2 (S,R 3 ) and n L 2 (S) , minimizing the worst case error means to minimize the sum of the two terms in the last inequality above. The solution to the extremal problem studied in the forthcoming sections offers a tool to trade off between them.

A bounded extremal problem

We consider the following bounded extremal problem (BEP for short):

Given e ∈ Ran b * 3 ⊂ L 2 (S, R 3 ) and a real number M > 0, BEP: find φ opt ∈ W 1,2 0 (Q) with ∇ φ opt L 2 (Q,R 2 ) ≤ M such that min φ∈W 1,2 0 (Q), ∇ φ L 2 (Q,R 2 ) ≤M b * 3 [φ] -e L 2 (S,R 3 ) = b * 3 [φ opt ] -e L 2 (S,R 3 ) .

Well posedness

Proposition 1 There exists a unique solution

φ opt to BEP; if e ∈ b * 3 W 1,2 0 (Q) , the constraint is saturated: ∇ φ opt L 2 (Q,R 2 ) = M . Proof: Since ∇φ L 2 (Q,R 2 ) is a norm equivalent to φ W 1,2 (Q) on W 1,2 0 (Q) by (3), the convex set {φ ∈ W 1,2 0 (Q) , ∇ φ L 2 (Q,R 2 ) ≤ M } is weakly compact in the Hilbert space W 1,2 0 (Q). Then, since b * 3 : W 1,2 0 (Q) → L 2 (S, R 3
) is continuous hence also weakly continuous [8, Thm 3.10], the set of approximants

C = b * 3 {φ ∈ W 1,2 0 (Q) , ∇ φ L 2 (Q,R 2 ) ≤ M } (16)
is weakly compact in L 2 (S, R 3 ). In particular it is weakly closed, and a fortiori it is closed in the norm topology. This implies there exists a unique best approximation to e from the closed convex set C, that can be written b

* 3 [φ opt ] for some φ opt ∈ W 1,2 0 (Q) which is unique because b *
3 is injective. This ensures both existence and uniqueness of the solution φ opt to BEP.

Next, assume that ∇ φ opt L 2 (Q,R 2 ) < M . In this case, the minimum value of the criterion is achieved at φ opt which lies interior to the approximation set. We then get by

differentiating b * 3 [φ] -e 2 L 2 (S,R 3 ) with respect to φ ∈ W 1,2 0 (Q) at φ opt that b * 3 [φ opt ] -e , b * 3 [δ φ ] L 2 (S,R 3 ) = b 3 b * 3 [φ opt ] -b 3 [e] , δ φ L 2 (Q) = 0 , ∀δ φ ∈ W 1,2 0 (Q) . ( 17
) Hence b 3 b * 3 [φ opt ] -b 3 [e] is orthogonal to W 1,2 0 (Q) in L 2 (Q), so by the density of W 1,2 0 (Q) in L 2 (Q) we must have that b 3 b * 3 [φ opt ] -b 3 [e] = 0
[φ opt ] ∈ b * 3 W 1,2 0 (Q) .
This contradicts the assumption on e, thereby achieving the proof. [START_REF] Jerison | The inhomogeneous Dirichlet problem in Lipschitz domains[END_REF] is not just closed but actually compact in L 2 (S, R 3 ). This furnishes another proof of the existence of φ opt , independent of convexity, which can be useful to deal with additional, possibly non-convex constraints in BEP. However, we do not consider such generalizations here.

Remark 1 The above proof shows that the constraint is saturated even if

e = b * 3 [φ] for some φ ∈ W 1,2 0 (Q), provided that ∇φ L 2 (Q,R 2 ) ≥ M . Remark 2 It follows from the compactness of b * 3 : W 1,2 0 (Q) → L 2 (S, R 3 ) that the set C defined in

The critical point equation

Let φ opt be the solution to BEP whose existence and uniqueness was proven in the previous section. Being the projection of e onto the closed convex set C ⊂ L 2 (S, R 3 ), the vector field b * 3 [φ opt ] is characterized by the variational inequality [8, Thm 5.2]:

e -b * 3 [φ opt ] , b * 3 [φ] -b * 3 [φ opt ] L 2 (S,R 3 ) ≤ 0, ∀φ ∈ W 1,2 0 (Q) , ∇φ L 2 (Q,R 2 ) ≤ M . ( 18
)
As we deal here with smooth constrained optimization, we can derive a more specific critical point equation (in short: CPE) to characterize φ opt . The CPE can be used to design numerical algorithms, and it also shows that φ opt is pointwise more regular than should a priori be expected from a W 1,2 0 (Q)-function, namely it is Hölder continuous of exponent 1/2. This substantiates a previous claim that constraining the derivative has a smoothing effect on our net moment estimator.

Proposition 2 Let e ∈ Ran b * 3 \ b * 3 W 1,2 0 (Q) ⊂ L 2 (S, R 3 ) and M > 0. i) A function φ opt ∈ W 1,2 0 (Q) is the solution to BEP iff ∇φ opt L 2 (Q,R 2 )
= M and there exists λ > 0 such that the following critical point equation holds, in the distributional sense on

Q: b 3 b * 3 [φ opt ] -λ Δ φ opt = b 3 [e] . ( 19
)
ii) The function φ opt lies in C ∞ (Q), and in the fractional Sobolev space W 3/2,2 (Q), as well as in C 1/2 (Q).

Proof:

If we let B : W 1,2 0 (Q) → L 2 (Q, R 2 ) be defined as B(φ) = ∇φ, then φ opt minimizes over φ ∈ W 1,2 0 (Q) the quantity b * 3 [φ] -e L 2 (S,R 3 ) subject to B(φ) L 2 (Q,R 2 ) ≤ M . Now, i) follows from [9, Thm 2.1] if we identify BB * (φ) with the distribution -Δφ, which is possible since C ∞ 0 (Q) is dense in W 1,2 0 (Q).
It is also instructive to establish [START_REF] Lang | Real and Functional Analysis[END_REF] directly, using differentiation. For this, let M ⊂ W 1,2 0 (Q) be the smooth hypersurface comprised of those φ such that ∇ φ L 2 (Q,R 2 ) = M . The tangent space T M φ to M at φ is the kernel of the linear form on W 1,2 0 (Q) given by ψ → ∇φ, ∇ψ L 2 (Q,R 2 ) . From Proposition 1 we know that φ opt ∈ M and, by optimality of φ opt , the derivative of φ → b * 3 [φ] -e 2 L 2 (S,R 3 ) at φ opt must vanish on T M φopt . Differentiating as we did in the proof of Proposition 1 and expressing that two linear forms with the same kernel must be proportional, we deduce there exists λ ∈ R such that, for every

ψ ∈ W 1,2 0 (Q), b * 3 [φ opt ] -e , b * 3 [ψ] L 2 (S,R 3 ) + λ ∇ φ opt , ∇ ψ L 2 (Q,R 2 ) = 0 . ( 20
)
Clearly we get an equivalent equation upon restricting ψ in [START_REF] Lima | Fast inversion of magnetic field maps of unidirectional planar geological magnetization[END_REF] to range over C ∞ 0 (Q), by density of the latter in W 1,2 0 (Q). Thus, by the Green formula for Sobolev functions, we find that ( 20) is equivalent to

b 3 b * 3 [φ opt ] -b 3 [e] , ψ L 2 (Q) -λ Δφ opt , ψ L 2 (Q) = 0 , ∀ψ ∈ C ∞ 0 (Q),
which means precisely that [START_REF] Lang | Real and Functional Analysis[END_REF] holds in the distributional sense. To see that λ > 0, substitute ψ = φ opt in [START_REF] Lima | Fast inversion of magnetic field maps of unidirectional planar geological magnetization[END_REF] to obtain:

b * 3 [φ opt ] -e , b * 3 [φ opt ] L 2 (S,R 3 ) = -λ ∇ φ opt 2 L 2 (Q,R 2 ) = -λ M 2 , ( 21
)
and observe from [START_REF] Kuttler | Eigenvalues of the Laplacian in two dimensions[END_REF], with φ = 0, that the above quantity is non-positive hence λ ≥ 0. Moreover λ = 0, otherwise [START_REF] Lima | Fast inversion of magnetic field maps of unidirectional planar geological magnetization[END_REF] would imply that ( 17) is satisfied and then, arguing as in the proof of Proposition 1, it would entail that e ∈ b * 3 W 1,2 0 (Q) , a contradiction. Conversely, assume that ∇φ opt L 2 (Q,R 2 ) = M and that [START_REF] Lang | Real and Functional Analysis[END_REF] holds for some λ > 0. Then Equation (20) holds as well, and so does Equation ( 21) which is a special case of [START_REF] Lima | Fast inversion of magnetic field maps of unidirectional planar geological magnetization[END_REF]. Subtracting ( 21) from [START_REF] Lima | Fast inversion of magnetic field maps of unidirectional planar geological magnetization[END_REF], we get that,

∀ψ ∈ W 1,2 0 (Q), b * 3 [φ opt ] -e , b * 3 [ψ] -b * 3 [φ opt ] L 2 (S,R 3 ) = -λ ∇ φ opt , ∇ ψ L 2 (Q,R 2 ) -M 2 ,
and by the Cauchy-Schwarz inequality we find that (18) holds when ∇ψ

L 2 (Q,R 2 ) ≤ M .
Hence φ opt is indeed the solution to BEP, thereby proving i).

As to point ii), let α = (α 1 , α 2 ) be a multi-index and set ∂ α to mean

∂ α 1 x 1 ∂ α 2 x 2 .
Since distributional derivatives commute, we see from [START_REF] Lang | Real and Functional Analysis[END_REF] that h := Δ ∂ α φ opt a fortiori belongs to L 2 (Q) because elements of the range of b 3 are restrictions to Q of real analytic functions on R 2 . Now, if we let p(z) = -Q log |z -t|h(t) dt be the logarithmic potential of χ Q h, we have that Δp = χ Q h and it is a standard consequence of the Calderòn-Zygmund theory that p ∈ W 2,2 (R 2 ) [2, Thm 4.5.3]. Then, ∂ α φ opt -p is harmonic hence real analytic in Q, implying that the restriction of ∂ α φ opt to any relatively compact disk D in Q lies in W 2,2 (D). Since α was arbitrary and W 2,2 (D) consists of continuous functions on D by the Sobolev embedding theorem [START_REF] Demengel | Espaces fonctionnels, Utilisation dans la résolution des équations aux dérivées partielles[END_REF]Ch. 4], we get that

φ opt ∈ C ∞ (Q).
Finally, as φ opt ∈ W 1,2 0 (Q) by definition, it has zero trace on the boundary of the bounded and Lipschitz-smooth domain Q, and since Δφ opt ∈ L 2 (Q) it follows from [16, Thm B, 2.] that φ opt belong to W 3/2,2 (Q). By the fractional version of the Sobolev embedding theorem, it follows that φ opt ∈ W 1,4 (Q) and consequently, by the standard Sobolev embedding theorem, it holds that φ opt ∈ C 1/2 (Q).

Remark 3

The restriction to exponent 3/2 made in point ii) on the Sobolev smoothness of φ opt is due to singularities that may occur on the Lipschitz boundary ∂Q. For instance if ∂Q is C ∞ -smooth, then by elliptic regularity theory we get that φ opt is C ∞ -smooth on Q, for Δφ opt is real analytic there and φ opt vanishes on ∂Q, see [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF].

We mention in passing that another application of [9, Thm 2.1] to the solution of extremal problems for harmonic gradients may be found in [START_REF] Atfeh | Bounded extremal and Cauchy-Laplace problems on the sphere and shell[END_REF]Prop. 4].

Analysis of the CPE and resolution schemes 4.1 Dependence on the constraint and the Lagrange parameter

The easiest way to make use of ( 19) is to pick λ > 0 and to solve for φ opt . This can be done in several manners, two of which will be discussed in this section. However, we no longer have control on the level of constraint M = ∇φ opt L 2 (Q,R 2 ) when doing this. This is why we need to know more about the behavior of φ opt and M as functions of λ. Hereafter, we use φ opt (λ) and M (λ) = ∇φ opt (λ) L 2 (Q,R 2 ) when we want to make the dependence on λ explicit.

Lemma 1 Let the assumptions be as in Proposition 2. In [START_REF] Lang | Real and Functional Analysis[END_REF], the following three statements are equivalent: i) the level of constraint M (λ) tends to +∞;

ii) the criterion b * 3 [φ opt (λ)] -e L 2 (S,R 3 ) tends to 0;

iii) the Lagrange parameter λ tends to 0.

Moreover, M (λ) is strictly decreasing as a function of λ, and M (λ) → 0 when λ → +∞.

Both φ opt (λ) and M (λ) are C ∞ -smooth functions of λ ∈ (0, +∞), with values in W 1,2 0 (Q) and R + respectively.

As an illustration, let us mention that the curves obtained numerically in Section 5 and shown in Figures 2 and3 provide a good overview of the properties listed in the lemma.

Proof: i) ⇒ ii): Assume that M → +∞. Since e ∈ b * 3 W 1,2 0 (Q) , the minimizing property of φ opt ensures that b * 3 [φ opt ] -e L 2 (S,R 3 ) → 0. ii) ⇒ iii): if b * 3 [φ opt ] -e L 2 (S,R 3 ) → 0, φ opt L 2 (Q) cannot become small since e = 0 and b * 3 is continuous. Now, observe that Δφ opt L 2 (Q) φ opt L 2 (Q) ≥ | Δφ opt , φ opt L 2 (Q) | = M 2 ≥ 1 C 2 φ opt 2 L 2 (Q)
by the Cauchy-Schwarz inequality, the Green formula and the Poincaré-Sobolev inequality. Therefore, Δφ opt L 2 (Q) cannot become small either. Because [START_REF] Lang | Real and Functional Analysis[END_REF] implies that

λ Δφ opt L 2 (Q) = b 3 [b * 3 [φ opt ] -e] L 2 (Q) ≤ b 3 b * 3 [φ opt ] -e L 2 (S,R 3 ) it necessarily holds that λ → 0. iii) ⇒ i): assume that λ → 0 but nevertheless M = ∇ φ opt L 2 (Q,R 2 ) → ∞.
Then, there exists a sequence λ k → 0 for which M (λ k ) is bounded, and we see from Equation ( 20)

that b * 3 [φ opt (λ k )] -e , b * 3 [ψ] L 2 (S,R 3 ) → 0 , ∀ψ ∈ W 1,2 0 (Q) . This implies that b * 3 [φ opt (λ k )] converges weakly to e in L 2 (S, R 3 ): indeed, consider Γ ∈ L 2 (S, R 3
) and an arbitrary ε > 0 ; we decompose Γ as Γ = Γ + Γ ⊥ where Γ ∈ D S and

Γ ⊥ ∈ D ⊥ S . Since b * 3 W 1,2 0 (Q) is dense in D ⊥ S = Ran b * 3 by (12), there exists ψ ∈ W 1,2 0 (Q) such that b * 3 [ψ] -Γ ⊥ L 2 (S,R 3 ) ≤ ε/ e L 2 (S,R 3
) and hence:

| b * 3 [φ opt (λ k )] -e, Γ L 2 (S,R 3 ) | ≤ | b * 3 [φ opt (λ k )] -e, Γ ⊥ -b * 3 [ψ] L 2 (S,R 3 ) | + | b * 3 [φ opt (λ k )] -e, Γ L 2 (S,R 3 ) | + | b * 3 [φ opt (λ k )] -e, b * 3 [ψ] L 2 (S,R 3 ) | ≤ ε + | b * 3 [φ opt (λ k )] -e, b * 3 [ψ] L 2 (S,R 3 ) |,
where we applied that b * 3 [φ opt ] -e L 2 (S,R 3 ) ≤ e L 2 (S,R 3 ) (because 0 is a candidate approximant) together with the Cauchy-Schwarz inequality to bound the first term, and we used the fact that b * 3 [φ opt ] -e belongs to D ⊥ S (cf. proof of Proposition 1), hence is orthogonal to Γ so that the second summand in the middle term of the inequality is zero. Finally, for k large enough (depending only on Γ and ε), the last term is smaller than ε, which proves the announced weak convergence. Now, rearranging [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF] as To prove that λ → φ opt (λ) is smooth, we introduce the operator T :

b * 3 [φ opt (λ k )] -e 2 L 2 (S,R 3 ) = -b * 3 [φ opt (λ k )] -e , e L 2 (S,R 3 ) -λ k M 2 (λ k ) , we find since λ k M 2 (λ k ) → 0
L 2 (Q) → W 1,2 0 (Q) that maps ψ ∈ L 2 (Q) to the solution u ∈ W 1,2 0 (Q) of the distributional equation Δu = ψ.
To see that T is well-defined, let us introduce the bilinear form a 0 defined for u, v

∈ W 1,2 0 (Q) by a 0 (u , v) = ∇u , ∇v L 2 (Q,R 2 ) . By density of W 1,2 0 (Q) in L 2 (Q), the equation Δu = ψ holds if and only if Δu , v L 2 (Q) = ψ , v L 2 (Q) for all v ∈ W 1,2 0 (Q),
which is equivalent (thanks to the Green formula) to a 0 (u , v) = -ψ , v L 2 (Q) . Now, the form a 0 is continuous and coercive by Equation ( 3), therefore the Lax-Milgram theorem (cf. [8, Cor. 5.8]) ensures the existence and uniqueness of a solution T (ψ) := u ∈ W 1,2 0 (Q). Notice, using successively the Cauchy-Schwarz inequality and Equation (2), that

∇u 2 L 2 (Q,R 2 ) = |a 0 (u, u)| = -ψ , u L 2 (Q) ≤ ψ L 2 (Q) u L 2 (Q) ≤ C ψ L 2 (Q) ∇u L 2 (Q,R 2 ) , whence ∇u L 2 (Q,R 2 ) ≤ C ψ L 2 (Q) . By Equation (3) we see that T (ψ) W 1,2 (Q) = u W 1,2 (Q) is bounded by C √ C 2 + 1 ψ L 2 (Q)
, which shows that T is indeed continuous from L 2 (Q) to W 1,2 0 (Q). More in fact is true, as hinted at in the proof of Proposition 2, namely T maps continuously

L 2 (Q) into W 3/2,2 0 (Q) ⊂ W 1,2
0 (Q). This is really the meaning of [16, Thm B, 2.], which says in our notation: "if ψ ∈ L 2 (Ω), then u ∈ W 3/2,2 (Ω)". Still, it seems difficult to reference an explicit statement, or even a formal argument. As explained in [16, p. 165], and also in the proof of Proposition 2, point ii), the logarithmic potential p of χ Q ψ satisfies Δp = ψ in Ω and meets an estimate of the form p W 2,2 (Ω) ≤ c ψ L 2 (Q) , but its trace on the boundary ∂Ω is not zero. Hence, u is obtained by subtracting from p a harmonic function in W 1,2 (Ω) whose trace is the same as the trace of p. The latter is easily checked to lie in W 1,2 (∂Ω) (defining Sobolev spaces on Lipschitz manifolds is possible because Lipschitz changes of variable preserve Euclidean Sobolev spaces, cf. [26, Thm 2.2.2]), therefore one is left to prove that a harmonic function u

∈ W 1,2 (Ω) with trace θ ∈ W 1,2 (∂Ω) satisfies u W 3/2,2 (Ω) ≤ c θ W 1,2 (∂Ω)
when Ω is Lipschitz smooth. This may be seen by combining the estimate [11, Eq. (1.1)], putting Φ(x) = x 2 therein and integrating on ∂Ω, with [16, Thm 4.1] which gives an equivalent quantity for the W 3/2,2 (Ω)-norm of a harmonic function. This last reference assumes the dimension is at least 3, but the planar case (which is our concern here) may be treated in the same way. The result we need is also explicitly stated and formally proved in [4, Prop. 5.2], which deals exclusively with dimension 2 (but more general equations). To sum up, T is not just continuous but also compact L 2 (Q) → W 1,2 0 (Q), by compactness of the embedding [START_REF] Demengel | Espaces fonctionnels, Utilisation dans la résolution des équations aux dérivées partielles[END_REF]Thm 4.54].

W 3/2,2 0 (Q) → W 1,2 0 (Q), see
We now turn to the proof that λ → φ opt (λ) is smooth. To this end, we first rewrite [START_REF] Lang | Real and Functional Analysis[END_REF] in the form F (λ, φ opt ) = 0 with F : R + × W 1,2 0 (Q) → W 1,2 0 (Q) given by

F (λ, φ) = φ - T (b 3 b * 3 [φ]) -T (b 3 [e]) λ , (22)
Differentiating ( 22) with respect to φ, which is possible since it is affine in φ, we get that D φ F (λ, φ) = I -T b 3 b * 3 /λ, where I denotes the identity operator. Now, T b

3 b * 3 : W 1,2 0 (Q) → W 1,2
0 (Q) is the composition of four operators, namely the embedding W 1,2 0 (Q) → L 2 (Q), the operator b * 3 , the operator b 3 and the operator T . Therefore T b 3 b * 3 is compact, for it is enough to reach the conclusion that one of these operators be compact, whereas all of them are! Being a compact perturbation of the identity, D φ F (λ, φ) will be invertible if only it is injective by a classical result of F. Riesz, see e.g. [

u ∈ W 1,2 0 (Q) satisfies T b 3 b * 3 [u]
= λu then we get from the Green formula, and since Δu ∈ L 2 (Q), that

-∇u 2 L 2 (Q,R 2 ) = Δu, u L 2 (Q) = λ -1 b 3 b * 3 [u] , u L 2 (Q) = λ -1 b * 3 [u] 2 L 2 (S,R 3 )
which implies ∇u = 0 since the first term is non-positive and the last one non-negative, hence also u = 0 by [START_REF] Astala | Elliptic Partial Differential Equations and Quasiconformal Mappings in the Plane[END_REF]. In particular D φ F (λ, φ opt (λ)) is an isomorphism of W 1,2 0 (Q) for each λ > 0 and so, by the implicit function theorem in Banach spaces [19, Ch. XIV, Thm 2.1], we deduce that λ

→ φ opt (λ) is C ∞ -smooth on R + since F : R + × W 1,2 0 (Q) → W 1,2 0 (Q) is obviously smooth. Specifically, φ opt (λ) satisfies a linear non-autonomous differential equation whose dynamics operator is the resolvent of T b 3 b * 3 : dφ opt dλ = (λI -T b 3 b * 3 ) -1 T (b 3 [e] -b 3 b * 3 [φ opt ]) λ = -(λI -T b 3 b * 3 ) -1 φ opt , ( 23
)
where the first equality holds because D λ F (λ, φ opt (λ)) + D φ F (λ, φ opt (λ))( dφopt(λ) dλ ) = 0 and the second uses [START_REF] Lang | Real and Functional Analysis[END_REF] and the fact that T Δφ opt = φ opt since Δφ opt ∈ L 2 (Q).

Next, we observe that for u ∈ W 1,2 0 (Q) with u = 0, the Green formula and the fact that ΔT = I on L 2 (Q) together imply:

∇(λI -T b 3 b * 3 )u , ∇u L 2 (Q,R 2 ) = λ ∇u 2 L 2 (Q,R 2 ) + ΔT b 3 b * 3 [u] , u L 2 (Q) = λ ∇u 2 L 2 (Q,R 2 ) + b 3 b * 3 [u] , u L 2 (Q) = λ ∇u 2 L 2 (Q,R 2 ) + b * 3 [u] , b * 3 [u] L 2 (Q,R 3 ) > 0, and setting u = (λI -T b 3 b * 3 ) -1 v gives us ∇(λI -T b 3 b * 3 ) -1 v , ∇v L 2 (Q,R 2 ) > 0, ∀v ∈ W 1,2 0 (Q), v = 0. ( 24
)
Now, the smoothness of the W 1,2 0 (Q)-valued function λ → φ opt (λ) entails that λ → ∇φ opt (λ) is continuously differentiable with values in L 2 (Q, R 2 ), and that d(∇φ opt )/dλ = ∇(dφ opt /dλ). Thus, it follows from ( 23) and ( 24) that

dM 2 dλ = d dλ ∇φ opt , ∇φ opt L 2 (Q,R 2 ) = 2 ∇ d dλ φ opt , ∇φ opt L 2 (Q,R 2 ) = -2 ∇(λI -T b 3 b * 3 ) -1 φ opt , ∇φ opt L 2 (Q,R 2 ) < 0,
thereby showing that M 2 (λ) is smooth and strictly decreasing as a function of λ. Obviously the same holds for M (λ), since it is never zero by the first assumption and the continuity of λ → M (λ). Finally, since M (λ) remains bounded as λ → +∞ by the first assertion of the lemma (or the monotonicity just proven), it follows from (2) that φ opt L 2 (Q) remains bounded as well. Then [START_REF] Lang | Real and Functional Analysis[END_REF] entails that Δφ opt → 0 in L 2 (Q) as λ → +∞. Applying the operator T now shows that φ opt (λ) → 0 in W 1,2 0 (Q) when λ → +∞, hence M (λ) → 0.

The monotonicity of λ → M (λ) asserted in Lemma 1 provides us with a means to approximately solve BEP if we can compute φ opt (λ), because we can use dichotomy on λ. Therefore we are reduced to estimate φ opt (λ) numerically for fixed λ. It is natural to ask if the fairly explicit differential equation in [START_REF] Stein | Singular Integrals and Differentiability Properties of Functions[END_REF] can be used for this purpose, however several issues arise in this connection. For instance the obvious boundary condition is φ opt (λ) → 0 in W 1,2 0 (Q) when λ → +∞, but the dynamics is singular at infinity. Moreover, integrating backwards with respect to λ is unstable. In addition, the equation is infinite-dimensional and a natural basis in which to truncate φ opt consists of eigenvectors of T b 3 b * 3 , which seem hardly accessible. In the forthcoming sections, we illustrate how the previous results can be used to derive a net moment estimator, but we rely on a direct approach based on numerical integration of (19).

Local moments

For S an open set in R 2 such that S ⊂ S, we define the local moment of the magnetization m on S to be m S = S m. For magnetizations with a density like those considered in this paper (i.e. m ∈ L 2 (S, R 3 ) in our case), recovering the local moment on any disk D(a, r) is tantamount to recovering m because by Lebesgue's monotone convergence theorem m D(a,r) /(πr 2 ) → m(a) for a.e. a ∈ S when r → 0. Hence it is clear from the outset that local moments cannot be recovered from the field in general, otherwise it would contradict the existence of nonzero silent magnetizations. Still the situation is more nuanced than it seems, and worth a short discussion.

First, since vector fields in L 2 (S, R 3 ) with vanishing horizontal component are orthogonal to D S , by [START_REF] Demengel | Espaces fonctionnels, Utilisation dans la résolution des équations aux dérivées partielles[END_REF], we may set e = (0, 0, χ S ) t ∈ D ⊥ S in BEP and the results we obtained so far apply to any S . Thus, the numerical procedures described in sections to come can be used in principle to estimate the local moments of the vertical component of m. This is consistent with the fact that compactly supported silent magnetizations have null vertical component.

The situation for the vector fields (χ S , 0, 0) t and (0, χ S , 0) t is different, for they do not belong to D ⊥ S . However, if V is an open neighborhood of S such that V ⊂ S, there exists a C 1 -smooth function ϕ on S which is identically 1 on S and 0 on S \ V . Then, the functions

ψ 1 (x 1 , x 2 ) = ϕ(x 1 , x 2 )x 1 and ψ 2 (x 1 , x 2 ) = ϕ(x 1 , x 2 )
x 2 are supported in V and lie in W 1,2 0 (S), so the gradients ∇ψ 1 , ∇ψ 2 do lie in D ⊥ S and they coincide respectively with (χ S , 0, 0) t and (0, χ S , 0) t on S . Now, if ∂S does not intersect the support of m, by compactness of the latter we can pick V to contain no points of supp m except those contained in S , and then S ∇ψ i . m = m i S for i = 1, 2. That is to say, if we know an open set S with S ⊂ S whose boundary does not intersect supp m (which means of course we have some a priori knowledge on m), we can construct an estimator of the local moment m S .

In another connection, whereas (χ S , 0, 0) t and (0, χ S , 0) t do not belong to D ⊥ S , their projections e 1,S and e 2,S onto that space in L 2 (S, R3 ) can be taken to be e in BEP. Then, by [START_REF] Demengel | Espaces fonctionnels, Utilisation dans la résolution des équations aux dérivées partielles[END_REF] and the definition of the orthogonal projection, solving this BEP will provide us with an estimator of the local moment on S of the magnetization of minimum L 2 (S, R 3 )norm which produces the same field as m. Note that e 1,S and e 2,S can be computed from the Hodge decomposition of (χ S , 0, 0) t and (0, χ S , 0) t by solving a Neumann problem, identical to the one set up in the proof of [START_REF] Baratchart | Silent and equivalent magnetic distributions on thin plates[END_REF]Prop. 5].

Local moment estimation is a natural extension of our method and has important practical applications. However, a detailed analysis on this topic is outside the scope of this paper and will be addressed in a future work.

Solving procedures

Equation ( 19) has a unique solution φ opt ∈ W 1,2 0 (Q) given e ∈ L 2 (S, R 3 ) and λ > 0. Indeed, taking the scalar product with an arbitrary ψ ∈ W 1,2 0 (Q), we may use the Green formula to obtain the equivalent equation:

a b (ψ, φ opt ) = ψ, b 3 [e] L 2 (Q) , ∀ψ ∈ W 1,2 0 (Q) , ( 25
)
where the bilinear form

a b on [W 1,2 0 (Q)] 2 is defined for φ , ψ ∈ W 1,2 0 (Q) by: a b (ψ , φ) = b * 3 [ψ] , b * 3 [φ] L 2 (S,R 3 ) + λ ∇ψ , ∇φ L 2 (Q,R 2 ) .
Observing that

λ ∇φ 2 L 2 (Q,R 2 ) ≤ b * 3 [φ] 2 L 2 (S,R 3 ) + λ ∇φ 2 L 2 (Q,R 2 ) ≤ b 2 φ 2 L 2 (Q,R 2 ) + λ ∇φ 2 L 2 (Q,R 2 ) ,
where b is the bound defined at Equation [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], we see that a b is continuous and coercive:

c b φ 2 W 1,2 0 (Q) ≤ a b (φ , φ) ≤ C b φ 2 W 1,2 0 (Q) ,
where we can take, e.g., c b := λ/(C 2 + 1) and C b := max(b 2 , λ), with C the constant of the Sobolev-Poincaré inequality associated to Q (cf. Equation ( 2)). Therefore, appealing to the Lax-Milgram theorem [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]Cor. 5.8], Equation ( 25) admits a unique solution φ opt ∈ W 1,2 0 (Q). Under the additional assumption that e ∈ Ran b *

3 \b * 3 W 1,2 0 (Q) ⊂ L 2 (S, R 3
), Proposition 2 tells us that φ opt coincides with the solution of BEP corresponding to the level of constraint

M = ∇ φ opt L 2 (Q,R 2 ) .
A standard approach to numerically solve equations like ( 25) is to estimate the solutions in weak form, in a space V of finite dimension. So, let (ψ i ) i∈I be a finite family of functions 3 spanning a subspace V = Span {ψ i , i ∈ I} of W 1,2 0 (Q). Instead of Equation ( 25), we consider the restricted equation

a b (ψ, φ) = ψ, b 3 [e] L 2 (Q) , ∀ψ ∈ V, ( 26
)
where the solution φ is searched in the space V : φ = j∈I α j ψ j for some real-valued coefficients (α j ) j∈I . Considering Equation [START_REF] Ziemer | Weakly Differentiable Functions. Sobolev Spaces and Functions of Bounded Variation[END_REF] with ψ = ψ i for every i ∈ I then yields

the system A α = ψ i , b 3 [e] L 2 (S,R 3 ) i∈I
where A = A(λ) is the matrix with coefficients

A i,j = a b (ψ i , ψ j ), i, j ∈ I. ( 27
)
Solving this system yields the coefficients of the solution φ ∈ V of Equation ( 26), which, by construction, is such that a b (φ -φ opt , ψ) = 0 for any ψ ∈ V . Letting φ indicate the orthogonal projection of φ opt onto V and φ ⊥ = φ opt -φ indicate its orthogonal complement, we thus get

|a b (φ opt -φ, φ opt -φ)| = |a b (φ opt -φ, φ ⊥ + (φ -φ))| = |a b (φ opt -φ, φ ⊥ )| ≤ a b (φ opt -φ, φ opt -φ)| 1/2 |a b (φ ⊥ , φ ⊥ )| 1/2 ,
where the inequality is obtained by Cauchy-Schwarz theorem applied to the positive form a b . This reduces to

|a b (φ opt -φ, φ opt -φ)| 1/2 ≤ |a b (φ ⊥ , φ ⊥ )| 1/2
, whence, with the coercivity of a b :

φ opt -φ W 1,2 0 (Q) ≤ C b c b φ ⊥ W 1,2 0 (Q) . Notice that all functions v in V satisfy φ ⊥ W 1,2 0 (Q) ≤ v -φ opt W 1,2 0 (Q)
by Pythagorean theorem. This irreducible error is inherent in the choice of the finite dimensional space V for representing an approximate solution. Of course, one needs to ensure a priori that φ ⊥ W 1,2 0 (Q) is small enough. Typically, the space V will be controlled by some parameter ε such that the distance of any given function to V = V ε tends to zero when ε → 0. In our examples of Section 5, for instance, ε corresponds to the step size of the chosen mesh grid on Q.

In practice, the overall process is as follows. We fix a value ε > 0, so that the dimension of V ε is large while keeping the computational burden acceptable. Since we do not know the value of the Lagrange parameter λ associated with our desired level of constraint M , we iterate the following procedure, starting with some initial value λ = λ 0 > 0: compute the coefficients α j = α j (λ) for j ∈ I, together with the corresponding approximation φ to φ opt , and the associated constraint level ∇ φ L 2 (Q,R 2 ) . If the latter is within a satisfactory range of M , we stop. Otherwise we bisect with respect to variable λ, according to the rule indicated by Lemma 1: increasing it when the constraint level is too high, and decreasing it otherwise.

If at some point monotonicity fails numerically, it means that the computations are inaccurate, e.g., that φ opt is still fairly far from V ε (i.e., φ ⊥ W 1,2 0 (Q) is fairly large) and we should try a smaller value of ε.

The above-described procedure is simple and yields fairly good results on the synthetic examples reported in Section 5. However, the precision strongly depends on φ ⊥ W 1,2 0 (Q) and the computational burden on dim(V ε ). It is therefore important to pick a family (ψ i ) able to approximate φ opt with as few elements as possible, in spite of its oscillatory behavior in certain regions that can be guessed from Lemma 1 and confirmed numerically. For the experiments presented in the next section to illustrate the technique, we did not try to optimize the design of (ψ i ) and did favor simplicity. More sophisticated choices of the basis are left here for further study.

The previous approach aims at approximating φ opt in W 1,2 0 (Q)-norm but tells us nothing about pointwise convergence, in particular it does not reflect that φ opt ∈ C 1/2 (Q), see Proposition 2. This is why it seems worth describing another algorithm to solve Equation [START_REF] Lang | Real and Functional Analysis[END_REF], which is more demanding computationally in general but offers some guarantee in this respect. It is based on successive approximations of the CPE itself by standard Dirichlet problems. More precisely, for e ∈ L 2 (S, R 3 ), λ > 0, > 0, φ 1 ∈ W 1,2 0 (Q) and n ≥ 2, we consider the sequence of equations:

(I -λ Δ) φ n = (I -b 3 b * 3 ) φ n-1 + b 3 [e] , (28)
where φ n ∈ W 1,2 0 (Q) is the unknown and φ n-1 acts as a parameter which was computed at the previous step.

Let us first show that, given φ n-1 ∈ W 1,2 0 (Q), there exists a unique solution

φ n ∈ W 1,2 0 (Q) to (28). Indeed, the bilinear form on [W 1,2 0 (Q)] 2 defined by a(φ , ψ) = φ , ψ L 2 (Q) + λ ∇φ , ∇ψ L 2 (Q,R 2 ) = (I -λΔ)φ , ψ L 2 (Q)
is continuous and coercive, while taking the scalar product with an arbitrary ψ ∈ W 1,2 0 (Q) on both sides of (28) yields the equivalent formulation:

a(φ n , ψ) = (I -b 3 b * 3 ) φ n-1 + b 3 [e] , ψ L 2 (Q) , ∀ψ ∈ W 1,2 0 (Q),
which has a unique solution φ n ∈ W 1,2 0 (Q) by the Lax-Milgram theorem. Thus, given some initial

φ 1 ∈ W 1,2 0 (Q), we can define inductively a sequence of W 1,2 0 (Q)-functions (φ n ) n∈N *
, where φ n is the solution to (28) for n ≥ 2.

Proposition 3 Let e ∈ Ran b * 3 \ b * 3 W 1,2 0 (Q) ⊂ L 2 (S, R 3
), λ > 0, and φ 1 ∈ W 1,2 0 (Q). i) For > 0 small enough, (28) defines a sequence of functions (φ n ) for n ≥ 2 that converges in W 1,2 0 (Q) to the unique solution φ opt (λ) of the critical point equation [START_REF] Lang | Real and Functional Analysis[END_REF].

In particular ∇ φ n L 2 (Q,R 2 ) → M (λ) as n → ∞. ii) Actually, φ n ∈ C 1/2 (Q) for n ≥ 2 and (φ n ) converges to φ opt in C 1/2 (Q).
Proof: Multiplying ( 19) by and subtracting it from (28), we obtain that

(I -λ Δ) [φ n -φ opt ] = (I -b 3 b * 3 ) [φ n-1 -φ opt ] , (29)
and taking the scalar product with φ n -φ opt in L 2 (Q) yields:

φ n -φ opt 2 L 2 (Q) + λ ∇ (φ n -φ opt ) 2 L 2 (Q,R 2 ) = (I -b 3 b * 3 ) [φ n-1 -φ opt ] , φ n -φ opt L 2 (Q) ≤ I -b 3 b * 3 φ n-1 -φ opt L 2 (Q) φ n -φ opt L 2 (Q) . ( 30
)
Clearly we may assume that φ n = φ opt for all n ≥ 1, for if φ n 0 = φ opt it follows from (30) that φ n = φ opt for all n ≥ n 0 and there is nothing to prove. Now, using Equation (2) on the left-hand side and dividing on both sides by φ n -φ opt we get

φ n -φ opt L 2 (Q) ≤ I -b 3 b * 3 1 + λ C 2 φ n-1 -φ opt L 2 (Q) . ( 31
)
to the computation of the right-hand side of Equation [START_REF] Ziemer | Weakly Differentiable Functions. Sobolev Spaces and Functions of Bounded Variation[END_REF]. However, thanks to the properties of T , these expansions do not only converge in L 2 (Q), but indeed in C 1/2 (Q), while the first approach approximates φ opt only weakly. A thorough study of the precision required in the computation to make the C 1/2 (Q) convergence effective, though, is beyond the scope of the present paper.

Numerical aspects and illustrations

We ran preliminary numerical experiments on the direct resolution scheme proposed in Section 4.3, using a family of finite elements in W 1,2 0 (Q). Observe that Equation ( 19) is an elliptic partial differential-integral equation to be solved on a square Q, subject to a homogeneous Dirichlet boundary condition. We thus make use of a family (ψ p,q ) of bilinear rectangular elements on a square mesh of Q. These elements are piecewise affine in the variables x 1 , x 2 separately, and are the simplest ones to expand W 1,2 0 (Q)-functions on a rectangle, see [10, Ch. 2] and [START_REF] Ern | Theory and Practice of Finite Elements[END_REF].

Details of our implementation

Mesh on Q. We choose a number P ∈ N that controls the number of elements that we use to mesh the rectangle Q. In order to simplify the presentation, we take Q to be the square [-R, R] × [-R, R] × {h} ⊂ R 2 × {h}, and we define δ = 2R/(P + 1) and κ i = -R + iδ for i = 1, . . . , P , so that the points of coordinates (κ p , κ q ), 1 ≤ p, q ≤ P constitute the P × P interior nodes of a square mesh of Q with step δ.

We denote with Q p,q the square centered at the node (κ p , κ q ) with side-length equal to 2 δ and edges parallel to the axes. Observe that Q = ∪ p,q=1,••• ,P Q p,q but the union is not disjoint as the squares overlap (a generic Q p,q shares interior points with its 8 closest neighbors, Q p ,q such that p ∈ {p -1, p, p + 1} and q ∈ {q -1, q, q + 1}). This is illustrated on Figure 1.

Elements. We take for p, q = 1, • • • , P the W 1,2 0 (Q)-functions ψ p,q defined on Q by

ψ p,q (x 1 , x 2 ) = 1 - |x 1 -κ p | δ 1 - |x 2 -κ q | δ χ Qp,q .
By construction, each ψ p,q is 0 on the boundary and outside Q p,q (it actually belongs to W 1,2 0 (Q p,q ) ⊂ W 1,2 0 (Q)) and ranges from 0 to 1 inside Q p,q , taking value 1 exactly at its center (κ p , κ q ). These properties make (ψ p,q ) an interpolating family of Lagrange type (the so-called Q 1 type). In particular, a function φ ∈ Span{ψ p,q } 1≤p,q≤P , φ = 1≤p,q≤P α p,q ψ p,q , always satisfies φ(κ p , κ q ) = α p,q for all 1 ≤ p, q ≤ P . Indeed, if we extend the notations by taking κ 0 = -R, κ P +1 = R and α 0,q = α P +1,q = α p,0 = α p,P +1 = 0 for any 0 ≤ p, q ≤ P + 1, then φ is the unique function whose restriction to any [κ p , κ p+1 ] × [κ q , κ q+1 ] (with 0 ≤ p, q ≤ P ) is a function of the form (x 1 , x 2 ) → a + bx 1 + cx 2 + dx 1 x 2 (with a, b, c, and d real values depending of course on p and q) and satisfying φ(κ p , κ q ) = α p,q for any 0 ≤ p, q ≤ P + 1.

Figure 1: The mesh on Q, here with P = 4. The points of coordinates (κ p , κ q ) (1 ≤ p, q ≤ P ) are represented by bullets. The elementary squares Q p,q overlap, as shown in the diagram.

Integrals on S.

For S, we also take a square, namely S = [-s, s] × [-s, s] × {0} ⊂ R 2 × {0}. The computations require evaluating a number of integrals on S: this is done using the trapezoidal rule applied on a uniform N × N grid for some parameter N ≥ 3. Specifically, using the notations ρ = 2s/(N -1) and σ i = -s + iρ for i = 0, . . . , N -1, the integral of a function f on S is approximated by

S f ρ 2 N -1 i=0 N -1 j=0 w i,j f (σ i , σ j ),
where w 0,0 = w 0,N -1 = w N -1,0 = w N -1,N -1 = 1/4, w 0,j = w N -1,j = w i,0 = w i,N -1 = 1/2 for 1 ≤ i, j ≤ N -2 and w i,j = 1 otherwise. Dot products on Q. We occasionally need to compute the dot product of L 2 (Q) between a function φ of Span{ψ p,q } 1≤p,q≤P and a smooth function g defined on Q. When it happens, we first approximate g with the elements (ψ p,q ), i.e.,

g g approx = 1≤k,l≤P g(κ k , κ l ) ψ k,l ,
and we then evaluate exactly φ, g approx L 2 (Q) = Q φ g approx . More precisely, if φ = 1≤p,q≤P α p,q ψ p,q , we have

φ, g approx L 2 (Q) = 1≤p,q≤P 1≤k,l≤P α p,q g(κ k , κ l ) ψ p,q , ψ k,l L 2 (Q) . ( 35
)
Hence, the computation boils down to evaluating a discrete sum. The terms ψ p,q , ψ k,l L 2 (Q) are analytically precomputed, once and for all, from the explicitly known expressions of the elements (ψ p,q ). Notice that such a product is non-zero only when the interiors of Q p,q and Q k,l overlap and that its value depends only on p -k and q -l. Hence there are actually only a few integrals to compute.

Computation of the matrix A. In order to construct the matrix from Equation (27), we need to compute the dot products ∇ψ p,q , ∇ψ k,l L 2 (Q,R 2 ) and the dot products b

* 3 [ψ p,q ] , b * 3 [ψ k,l ] L 2 (S,R 3
) . Regarding the former, the gradients ∇ψ p,q are explicitly known from the above expressions for ψ p,q , so that these products can be analytically computed. The remarks we made for the computation of ψ p,q , ψ k,l L 2 (Q) apply here also.

As to the products b

* 3 [ψ p,q ] , b * 3 [ψ k,l ] L 2 (S,R 3
) they are computed as the sum of three integrals on S, each one being evaluated by the trapezoidal rule as explained above. This requires to numerically estimate b * 3 [ψ p,q ] (σ i , σ j ) for all 1 ≤ p, q ≤ P and 0 ≤ i, j ≤ N -1. We do it using Formula (9) with φ = ψ p,q : denoting with K the convolution kernel (K = ∂ x 1 P h or K = ∂ x 2 P h or K = (∂ x 3 P x 3 ) |x 3 =h depending on the component being considered), the value of (K φ)(σ i , σ j ) is evaluated, with the method described above, as the dot product

[K φ](σ i , σ j ) = K i,j , ψ p,q L 2 (Q) ,
where

K i,j : (x 1 , x 2 ) → K(σ i -x 1 , σ j -x 2 )
. The numerical evaluations of K i,j are obtained by an explicit formula for K i,j analytically derived from its definition and Equation [START_REF] Baratchart | Asymptotic method for estimating magnetic moments from field measurements on a planar grid[END_REF].

Right-hand side of the system. The right-hand side of Equation ( 26) requires to evaluate ψ p,q , b 3 [e] L 2 (Q) for all 1 ≤ p, q ≤ P . These are dot products on Q of the form described above, so we simply need to compute b 3 [e](κ k , κ l ) for all 1 ≤ k, l ≤ P . This is done using Equation ( 8): convolutions with the same kernels as before appear, except that the integrals involved are now on S and not on Q. We evaluate them with the trapezoidal rule described above.

Numerical results

Numerical simulations were performed using Matlab. The geometry of the measurement setup is fixed to R = 0.00255 and h = 0.00027, while the sample size is s = 0.00197. This corresponds to a small, yet realistic, sample studied by SQUID microscopy. The discretization parameters P and N have been chosen to be P = N = 100. We estimated the solution of the Critical Point Equation ( 19) using the above-described algorithm, for different values of λ, and for the right-hand sides corresponding to e = e k with k = 1 and k = 3. We denote with φ e k (λ) the approximation of the true solution φ opt so obtained. Notice that the case k = 2 is similar to the case k = 1, by symmetry (because φ e 2 (x 1 , x 2 ) = φ e 1 (x 2 , x 1 )), which is why we only show the results corresponding to e 1 and e 3 . Figures 2 and3 show how the criterion b * 3 [φ e k (λ)] -e k L 2 (S,R 3 ) and the level of constraint M (λ) = ∇φ e k (λ) L 2 (Q,R 2 ) depend on λ. This is in accordance with the assertions of Lemma 1 about the behavior of φ opt .

Figure 4 shows the same information in another form by plotting the evolution of the criterion with respect to the level of constraint. This is the so-called L-curve (see [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems[END_REF]): when the level of constraint M is fairly small, the criterion decreases very fast, meaning

kb $ 3 [? e1 (6)] ! e 1 k L 2 (S;R 3 ) ke 1 k L 2 (S;R 3 ) kb $ 3 [? e3 (6)] ! e 3 k L 2 (S;R 3 ) ke 3 k L 2 (S;R 3 ) Figure 2: Approximation error of b * 3 [φ e k (λ)
] with respect to e k when λ varies (the scale for λ is logarithmic). The solid line corresponds to the case when k = 1, while the dashed line corresponds to the case when k = 3. As expected from Lemma 1, on the one hand, when λ goes to 0 (i.e., log 10 (λ) → -∞), the error tends to 0. On the other hand, when λ goes large, the constraint M (λ) goes to 0, meaning that φ opt is forced to go to 0, whence the relative error tends to 1. that the quality of the linear estimator can be much improved, at the small cost of making it behave only slightly worse (a small increase of M intuitively means slightly larger oscillations of φ opt ). On the other hand, for larger values of M , the criterion is barely improved even when M is rather well increased, meaning that the benefit on the quality of the linear estimator is not worth the deterioration of its behavior. A compromise between both situation lies at the "elbow" of the L-curve, e.g., around values of M corresponding to λ = 10 -21 or λ = 10 -22 .

* 3 [φ e k (λ)] -e k L 2 (S,R 3 ) / e k L 2 (S,R 3 ) as a function of the constraint M = ∇φ e k (λ) L 2 (Q,R 2 )
. The upper plot corresponds to the case when k = 1, and the lower one to the case when k = 3. As expected, the error decreases and tends to 0 as the constraint is relaxed.

We show on Figure 5 the functions φ e k (λ) (k = 1 and 3) for λ = 10 -21 . They are plotted on a rectangle slightly larger than Q to help understand how they behave at the boundary of Q (they are, by definition, equal to 0 outside Q). On the bottom layer of each plot, a map is displayed where the color of each pixel corresponds to the value of the function. This does not say more than the plot itself, but in a complementary way, in order to help understanding the oscillations pattern of φ e k . Also, on the same bottom layer, the rectangles Q and S are displayed, so as to recall their respective positions. We can see that the functions φ e k already have a fairly important oscillation pattern when λ = 10 -21 , more specifically on the region corresponding to Q \ S. Indeed the maximal absolute value of φ e 1 on Q \ S is approximately 6.8 • 10 5 while its maximal absolute value on S is roughly 0.89 • 10 5 . The same observation holds for φ e 3 (with respective values 1.95 • 10 5 and 0.38 • 10 5 ). Also, we can see that the highest oscillations are close to the boundary of Q, and indeed the functions go to 0 fairly abruptly on ∂Q. This means that, when estimating the moment m from the data b 3 [m] with the linear estimator, much importance will be given to the values of the field measured near the edges, in spite of the fact that the vanishing of φ opt on ∂Q was imposed in hope to avoid this issue. The latter may be a drawback in practice, as measurements taken near the edges usually have lower signal-to-noise ratios, and further analysis is needed here to check how much this can affect the results with real data and whether there is a need to devise a remedy.

Finally, a close inspection of the plots, especially at some peaks of the functions, seems to suggest that the functions exhibit kinks and are not very smooth. This is likely to be an artifact due to the approximation of φ opt by a function in Span{ψ p,q } 1≤p,q≤P and indicates that the choice of P might be a bit too small to render the behavior of φ opt accurately at those places.

In order to test the ability of the computed estimators φ e k (λ) with λ = 10 -21 to recover net moments, we considered a synthetic magnetization m = (m 1 , m 2 , m 3 ) on S. The considered magnetization is shown on Figure 6 and is based on the logo of the Apics project-team. We discretized the logo on a 540 × 540 grid of dipoles. Each part of the logo (letter A, letters PICS, and mountain) is magnetized along a different direction. The dipoles belonging to each part have close moments, but not exactly equal, hence simulating an almost uniformly magnetized shape. Overall, the net moment of the 'A' part is approximately (-12, -86 , 3.5) • 10 -6 , the net moment of the 'PICS' part is approximately (-61, -26, 25) • 10 -6 and the net moment of the mountain is approximately (-0.76, -0.28, 13) • 10 -6 . The total net moment m of the synthetic magnetization is approximately (-74, -112, 41) • 10 -6 . The choice of using dipoles does not strictly match our framework of working with L 2 magnetizations, but the grid is intended to be fine enough that it might actually be viewed as a practical approximation of a piecewise continuous magnetization. The reason for this choice is practical: it allows us to simply use the exact formula for the field of a magnetic dipole in order to evaluate the forward operator b 3 .

We computed the values of b 3 [m] at the points of our mesh on Q. To test the influence of additive noise, we also generated a Gaussian white noise component with a magnitude of order roughly 1% of the maximal absolute value of b 3 [m]. Using these values, we approximated the field and the noise as functions on Q using the elements (ψ p,q ), in order to evaluate the dot products with functions φ e k (λ) (k = 1, 2, 3) as described in Equation (35). Figure 7 shows the corresponding functions, and Table 1 sums up our results.

In absence of noise, the relative error between the individual components m k and their estimates decrease as λ tends to 0 because b 3 [m] is smooth enough that its integral 6 and an additive Gaussian white noise generated on the same P × P mesh on Q, with P = 100. The computed values are used to approximate the field and the noise as functions of Span{ψ p,q } 1≤p,q≤P . Note that the color scales are not the same on both pictures.

against the linear estimators φ e k (λ) is fairly well evaluated. Interestingly, the component m 3 is more accurately recovered than the others, which is likely connected with the fact that the inverse problem of recovering m 3 from b 3 [m] has a unique solution (because silent sources have a null vertical component) whereas recovering m 1 or m 2 does not. This difference is reflected in the expression of the kernels generating ∂ x 1 Λ and ∂ x 2 Λ from m, compare ( 5) and [START_REF] Baratchart | Characterizing kernels of operators related to thin-plate magnetizations via generalizations of Hodge decompositions[END_REF].

In contrast, in the presence of noise, we see that the quality of the estimate is first getting better as λ decreases, but it reaches an optimum and the error starts increasing when λ goes below some level because the linear estimators peak so much and oscillate so fast that their dot product against the noise become non-trivial.

No noise

With noise λ δ 1 (%) δ 2 (%) δ 3 (%) δ r (%) θ ( • ) δ r (%) θ (

Concluding remarks

In this paper, we considered the inverse moment problem in magnetostatics, and derived linear estimators for the moment from the field in a planar geometry, We did this in a Hilbertian framework for L 2 magnetizations on thin slabs, and we studied a regularization technique based on the solution of a bounded extremal problem for the estimator. Encouraging numerical experiments, reported in Section 5, were carried out on a synthetic though already nontrivial example, with 1% synthetic Gaussian white noise added, using the simple resolution scheme described in Section 4.3. The treatment of real data by this method can still be enhanced by using more precise estimates of b * 3 [ψ p,q ] than those obtained using quadrature formulas, in order to get a more accurate version of the matrix A and of the right-hand side in [START_REF] Ziemer | Weakly Differentiable Functions. Sobolev Spaces and Functions of Bounded Variation[END_REF]. Moreover, finer (possibly adaptive) meshes can be used, though this will increase the computational burden. An alternating strategy, when the measurement area is a rectangle, rests on solving Equation (34) in the eigenbasis of the Laplacian, cf. the discussion at the end of Section 4.3.

However, when dealing with physical data, we should face additional problems arising from increased, non purely stochastic noise and other imperfections in the measurements. We did not touch upon this issue in the present paper whose main focus is theoretical, with numerical experiments designed for demonstrating the method under moderately realistic conditions.

It would be interesting to carry the approach over to magnetizations modeled by more general measures than those having L 2 -density, and to more regular (e.g., Lipschitzsmooth) estimators. This is needed to handle some popular models for magnetizations such as sums of dipoles. Although the problem looks more difficult, we expect that the method can be adapted to this case.

  cannot vanish for all m unless e i = b * 3 (φ i ), which is impossible because e i ∈ Ran b * 3 by [5, Lem. 7, (iii)]. On the one hand, unboundedness of µ entails that the moment recovery problem knowing b 3 [m] is ill-posed, in that small variations of b 3 [m] in the L 2 (Q)-metric may result in large variations of m . On the other hand, the fact that e i ∈ D ⊥ S = Ran b * 3

  . Thus, b * 3 [φ opt ] -e belongs to D S = Ker b 3 . However, both b * 3 [φ opt ] and e belong to D ⊥ S , see (12), thus b * 3 [φ opt ] -e = 0 and consequently e = b * 3

  by assumption and b * 3 [φ opt (λ k )] -e , e L 2 (S,R 3 ) → 0 by the weak convergence just proven that b * 3 [φ opt (λ k )] -e L 2 (S,R 3 ) → 0. But e ∈ b * 3 W 1,2 0 (Q) by hypothesis, therefore this limiting relation implies that M (λ k ) → +∞, as indicated at the end of Section 2.3. This contradiction achieves the proof of the first assertion.

Figure 3 :

 3 Figure3: Constraint M (λ) = ∇φ e k (λ) as a function of λ (the scale for λ is logarithmic). The solid line corresponds to the case when k = 1, while the dashed line corresponds to the case when k = 3. As expected from Lemma 1, these are strictly decreasing smooth functions tending to +∞ when λ → 0 (i.e., log 10 (λ) → -∞) and tending to 0 when λ → +∞.

20 6 = 10 ! 19 6 = 10 ! 18 6 = 10 ! 17 M 17 MFigure 4 :

 20101910181017174 Figure 4: "L-curves" showing the approximation error b* 3 [φ e k (λ)] -e k L 2 (S,R 3 ) / e k L 2 (S,R 3 ) as a function of the constraint M = ∇φ e k (λ) L 2 (Q,R 2 ). The upper plot corresponds to the case when k = 1, and the lower one to the case when k = 3. As expected, the error decreases and tends to 0 as the constraint is relaxed.

Figure 5 :

 5 Figure 5: φ e 1 (λ) (top) and φ e 3 (λ) (bottom) for λ = 10 -21 . On both plots, the rectangles Q (red) and S (blue) are drawn together on the bottom layer to help visualize their respective positions.

Figure 6 :

 6 Figure 6: Synthetic magnetization (from left to right) m 1 , m 2 and m 3 on S.

Figure 7 :

 7 Figure 7: Field b 3 [m] corresponding to the synthetic magnetization shown on Figure6and an additive Gaussian white noise generated on the same P × P mesh on Q, with P = 100. The computed values are used to approximate the field and the noise as functions of Span{ψ p,q } 1≤p,q≤P . Note that the color scales are not the same on both pictures.

  (k = 1, 2, 3) of the net moment m are approximated thanks to the linear estimator as µ k = b, φ e k (λ) L 2 (Q) for several values of λ and with b being either the exact synthetic field b 3 [m] or the exact field plus some noise. The quantities δ k are the relative errors δ k = µ k -m k / m k . The quantity δ r is the relative error of the amplitude of µ = (µ 1 , µ 2 , µ 3 ) as a vector approximation of m , i.e., δ r = µm / m , where • denotes the Euclidean norm. Finally, θ is the angle between the vectors µ and m , i.e., θ = 360 2π arccos µ µ . m m .

  of divergence free vector fields in L 2 (R 2 , R 2 ) by the extension theorem for Sobolev functions on a Lipschitz open set [12, Prop. 2.70].

	1,2 0 (Q)

⊥

are restrictions to S

  25, Ch. XVII, Prop. 2.3]. It is indeed injective at every argument (λ, φ) (note that it does not depend on φ), for if

Table 1 :

 1 The components m k

	• )

Mathematically speaking, m is the effect of the distribution m on any smooth compactly supported function on R 3 which is identically 1 in a neighborhood of S, see Section

2.2.

For instance, in our examples in Section 5 where we deal with square-shaped Q, we consider a mesh on Q with P ×P points, where P is a parameter. The set I then corresponds to {(p, q) ∈ N 2 , 1 ≤ p, q ≤ P }, and the functions (ψ i ) are Q 1 finite elements (see Section 5 for details).

Acknowledgments. All authors were supported in part by an Inria grant to the associate team Impinge. L. Baratchart, S. Chevillard, J. Leblond and E. A. Lima also acknowledge support from the MIT-France seed fund. The research of D. Hardin and the research of E. A. Lima were supported, in part, by the U.S. National Science Foundation under the grants DMS-1521749 and DMS-1521765, respectively. We are grateful to M. Northington and D. Ponomarev for many valuable discussions.

Next, observe that (b being the bound introduced with Equation ( 10))

Therefore, whenever 0 < < 1/b 2 , the operator I -b 3 b * 3 is positive self-adjoint on L 2 (Q), so its norm is (see e.g. [8, Prop. 6.9]):

which, in view of the equality in (32), is smaller than 1. From (31) we finally obtain

which establishes that φ n -φ opt L 2 (Q) decreases to 0 geometrically fast as n → ∞. Rewriting now (29) in the form:

where the operator T was defined in the proof of Lemma 1, we get from [16, Thm B, 2.] and the convergence of

Still, computing the solution to (28) with good accuracy on pointwise values is not an easy task. Let us simply mention that Equation (28) can be rewritten as the equivalent equation

An interesting feature of this formulation is that the factor T has a smoothing effect, converting L 2 (Q)-convergence into C 1/2 (Q)-convergence. Moreover, if Q is a rectangle, the eigenfunctions and eigenvalues of the Laplacian are explicitly known (the eigenfunctions are products of sines) [START_REF] Kuttler | Eigenvalues of the Laplacian in two dimensions[END_REF]. They form a Hilbert basis of L 2 (Q) on which the operators T and ( λI -T ) -1 are diagonal. Let us denote with (u p,q ) p,q∈N the eigenfunctions and assume that the expansion of φ n-1 = p,q∈N β p,q u p,q is known. Thanks to the orthonormality of the basis, the coefficients β k,l of the expansion of (

] in the basis are given by

Therefore, once the products b

Q) have been precomputed, it is easy to get the expansion of φ n from the expansion of φ n-1 . The computational burden in the case of a rectangle is overall fairly similar to the first approach, the precomputation of b * 3 [u p,q ], b * 3 [u k,l ] L 2 (S,R 3 ) playing the same role as the precomputation of the matrix A of Equation ( 27), and the precomputation of u p,q , b 3 [e] L 2 (Q) corresponding