MAGNETIC MOMENT ESTIMATION AND BOUNDED EXTREMAL PROBLEMS

. 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 𝐿 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.

1. 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 [5,6,7,20] 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 0 ⊂ R 3 , and is assumed to support a R 3 -valued magnetization distribution . The measurements consist of pointwise values, on a compact set Q contained in a plane ℎ parallel to 0 lying at distance ℎ > 0 from it, of the normal component of the magnetic field produced by (normal to ℎ , 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 × {ℎ}, respectively. Thus, the measured quantity is the vertical component 3 [ ] of the magnetic field produced by , evaluated on Q. Maxwell's equations for the magnetostatic case [15,Ch. 5] imply that this magnetic field derives from a scalar magnetic potential Λ which is a solution to the Poisson equation ΔΛ = div in R 3 . Consequently, we have that 3 [ ] = − 0 3 Λ with 0 = 4 × 10 −7 , see e.g. [7,5]. 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 [7], 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 [5] 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 3 [ ] = 0 on Q if and only if is a silent magnetization of the type described in [7]. 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 ⟨ ⟩ of the unknown magnetization , i.e. the integral 1 of . So, we merely aim at recovering a vector in R 3 (namely ⟨ ⟩) rather than the whole distribution on the rectangle S. Unlike magnetization recovery, net moment recovery is meaningful without extra-assumptions on , 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 [7,20] for more explanation.
From a physical viewpoint, knowing the net moment of a magnetization is a basic piece of information that is relevant in itself. It gains further significance when putting additional assumptions on , 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 [20] and [7,Sec. 4.1]. More can be said in this connection: in fact, to any magnetization compactly supported in a plane 0 , there is a unidimensional magnetization u supported on 0 which generates the same field on a given connected component of the complement of 0 . Surprisingly perhaps, the direction of u can even be picked arbitrarily provided it is not tangent to the plane [7,Thm 3.6]; however, clearly, u will not have the same support as in general, and will typically be non-compactly supported. Now, even in this case, choosing u to have the direction of ⟨ ⟩ 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 [7, Eqns (29)-(30))]. This dangles the prospect of efficiently computing a unidimensional magnetization which differs from the true one by a (generally noncompactly 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 3 [ ] are taken merely on Q which does not surround S. If measurements were available on the entire plane R 2 × {ℎ} (which does surround S if one takes into account points at infinity), then asymptotic formulas from [6] would allow us in principle to approximate the net moment arbitrary well, at least if 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 3 [ ] on Q may result in large differences on ⟨ ⟩, 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 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 for all ∈ 2 (S, R 3 ) of given norm, with a bound on the norm of the derivative of in 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 3 [ ] for estimating ⟨ ⟩, 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 3 [ ] 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 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 (1). Computing such local moments could be of interest to get information on the support of , 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.
2.1. Function spaces. For = ( 1 , · · · , ) ∈ R , we let ‖ ‖ = (︀ 2 1 + · · · + 2 )︀ 1 2 be its Euclidean norm; here and below, a superscript " " means "transpose". Denote the scalar product of and by . = ∑︀ . The notation stands irrespective of but no confusion should arise, and we are only concerned with = 2, 3. We write for the closure of a set and for = 2 we indicate the open disk of center and radius with ( , ). We write ∞ (Ω) for the space of smooth functions on an open set Ω ⊂ R having continuous derivatives of any order, and ∞ 0 (Ω) for the space of smooth functions with compact support in Ω. Recall that distributions on Ω are linear forms on ∞ 0 (Ω) which are continuous for a certain topology, the precise definition of which may be found in [22,Sec. I.2]. The support of a distribution d on Ω, denoted by supp d, is the largest closed set ⊂ Ω such that d( ) = 0 for all ∈ ∞ 0 (Ω ∖ ). Distributions get differentiated by the usual rule: d( ) = −d( ), with to mean the partial derivative with respect to the coordinate . It is standard to regard a locally integrable function as a distribution by putting ( ) = ∫︀ . For a distribution or a function, we let ∇ = ( 1 , · · · , ) for the gradient of . The divergence operator, which is the distributional adjoint to ∇, is denoted as ∇. and acts on R -valued distributions as ∇.( 1 , · · · , ) = Σ 1≤ ≤ .The Laplace operator Δ is defined by Δ = Σ 1≤ ≤ 2 . These operator notations stand irrespective of , but no confusion should arise. In fact we use them for = 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 2 (Ω, R ) the familiar Lebesgue space of R -valued square summable functions which is a Hilbert space with norm and scalar product: with dℓ to indicate Lebesgue measure. We simply write 2 (Ω) if = 1.
If Ω is open, we put 1,2 (Ω) for the standard Sobolev space of functions belonging to 2 (Ω) together with their first distributional derivatives. It is a Hilbert space with norm: We let If Ω is bounded, the Sobolev-Poincaré inequality [8,Cor. 9.19] implies there exists a constant > 0 (depending on Ω) such that Therefore ‖·‖ 1,2 (Ω) and ‖∇ ·‖ 2 (Ω,R 2 ) are equivalent norms on 1,2 0 (Ω), namely: When Ω is bounded and Lipschitz-smooth, meaning that its boundary Ω is locally the graph of a Lipschitz function, then functions in 1,2 (Ω) have a trace on Ω which belongs to the fractional Sobolev space 1/2,2 ( Ω), consisting of square summable functions with respect to arclength on Ω for which where |d | indicates (the differential of) arclength. The trace operator is continuous and surjective onto 1/2,2 ( Ω). Functions with zero trace are exactly those belonging to 1,2 0 (Ω), and they are also those whose extension by zero outside Ω defines a function in 1,2 (R 2 ). We also make use of the fractional Sobolev space 3/2,2 (Ω). The latter consists of functions ∈ 1,2 (Ω), each partial derivative of which satisfies: The Sobolev space 2,2 (Ω), consists of functions in 2 (Ω) whose partial derivatives of the first order lie in 1,2 (Ω). It is a Hilbert space with norm We refer to [1], [8,Ch. 9], [23, Ch. V, VI] for classical properties of Sobolev spaces. For 0 ≤ < 1, we put (Ω) for the space of Hölder continuous functions of exponent on Ω.

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 where the brackets denote the duality pairing between distributions and ∞ 0 (R 3 )functions. Note that the formula makes sense because → 1/‖ − ‖ coincides with a ∞ 0 (R 3 )-function on a neighborhood of supp m when / ∈ 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 = ( 1 , 2 , 3 ) there, then we get that Let 3 be the Poisson kernel of the upper half-space { 3 > 0}, see e.g. [24, Ch. I]: On this half-space, the previous expression for Λ implies that it is half the convolution over R 2 of 3 with the function 1̃︀ 1 + 2̃︀ 2 +̃︁ 3 , where 1 , 2 are the Riesz transforms on R 2 , see [7,Sec. 2]). In other words, for 3 > 0, 2Λ( 3 ) is the harmonic extension of 1̃︀ 1 + 2̃︀ 2 +̃︁ 3 to the upper half space. Likewise, for 3 < 0, 2Λ( 3 ) is the harmonic extension of 1̃︀ 1 + 2̃︀ 2 −̃︁ 3 to the lower half space.
Hereafter, we let and be two Lipschitz-smooth bounded non-empty open sets in R 2 × {0} and R 2 × {ℎ} respectively. For simplicity, we also assume that is simply connected. One may think of , as being rectangles whose closures , are the sets S, Q introduced in Section 1. Observe that expression (5) for Λ may as well be computed as an integral over because the boundary = S ∖ has Lebesgue measure zero. We often identify and with subsets of R 2 , while the third component (0 in the case of and ℎ in the case of ) is treated as a parameter.
The main operator under consideration here is to the vertical component of the magnetic field on . Specifically, 3 Although the target space is here 2 ( ), it is clear since ℎ > 0 that the range of 3 consists of restrictions to of smooth (even analytic) functions on { 3 > 0}. Basic properties of the operator 3 were given in [5,Sec. 3], [7]. Below, we recall some of them that will be used in the course of the paper.
Recall the notatioñ︀ = ∨ 0 ∈ 2 (R 2 ). From the relation we get for all It is easily checked that 3 is a bounded operator. In fact, being an integral operator whose kernel is uniformly continuous on × , it is compact [17, Ch. III, Ex. 4.1]. Its adjoint * 3 : 2 ( ) → 2 ( , R 3 ) acts on ∈ 2 ( ) by the formula and is likewise bounded, even compact because so is 3  ∨ 0 is silent. This means that 3 = 0 and that ( 1 , 2 ) belongs to the space ⊂ 2 ( , R 2 ) consisting of vector fields whose extension by zero outside of defines a divergencefree vector field on R 2 . In the proof of [5,Prop. 2], it is shown that can be parameterized by Sobolev functions: Besides, it follows at once from (11) and the density of ∞ 0 (Ω) in 1,2 0 (Ω) that the orthogonal space ⊥ to in 2 ( , R 2 ) is comprised of those vector fields satisfying the distributional Schwarz rule. Since is simply connected, these are the distributional gradients on , by [22, Ch. II, Sec. 6, Thm VI]. Thus, appealing to [12,Thm 6.74] which implies that a distribution with 2 derivative on a Lipschitz open set must be a Sobolev function 2 , we conclude that ⊥ is the set ∇ 1,2 ( ) of gradients of Sobolev functions. Altogether, this leads us to a characterization of the kernel of 3 and of the range of * 3 in 2 ( , R 3 ) as follows (compare [5,Lem. 4]). If we set = Ker 3 , then where ⊥ stands for the orthogonal space to in 2 ( , R 3 ) and the second equality in the second line of (12) holds because 1,2 0 ( ) is dense in 2 ( ). Also, rotating pointwise by /2 the relation ⊥ = ∇ 1,2 ( ), we get that All magnetizations ′ such that 3 [ ′ ] = 3 [ ] have the same net moment as . Indeed, 3 [ ] = 0 ⇔ ∈ , and since 1 , 2 are gradients of Sobolev functions on (namely of ↦ → 1 , ↦ → 2 ), we see from (12) that each element of is orthogonal to in 2 ( , R 3 ) for 1 ≤ ≤ 3. Hence, ⟨ ⟩ is uniquely determined by 3 [ ], in other words there is a linear map : . This map, however, is not continuous. Otherwise indeed, it would have a continuous extension 2 ( ) → R 3 by the Hahn-Banach theorem and thus, to each ∈ {1, 2, 3}, there would exist ∈ 2 ( ) such that the quantity vanishes for all ∈ 2 ( , R 3 ). But the last term in (13) Then, computing as in (13), we get that . Equation (14) shows we can estimate ⟨ ⟩ with arbitrary relative precision, at least in principle, by taking the scalar product of the data 3 [ ] with some appropriate estimator , . However we necessarily have that ‖ , ‖ 2 ( ) → +∞ when → 0, otherwise taking a weakly convergent subsequence would imply in the limit that ∈ Ran * 3 , a contradiction. Thus, when making small in order to improve accuracy of the estimate ⟨ 3 [ ], , ⟩ 2 ( ) for ⟨ ⟩, we tend to increase the norm of the estimator , 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, , we will control the 1,2 0 ( )-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 * 3 [ ]. Second, designing the estimator so as to vanish on the boundary puts smaller weight on the values of the field close to , which have poor signal/noise ratio because of its rapid decay when the distance to the sample increases. Last but not least, constraining the 1,2 0 ( )-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 ∈ Ran * 3 ⊂ 2 ( , R 3 ). In connection with net moment estimation, can be any of the for ∈ {1, 2, 3} and then we are in the generic case where ̸ ∈ Ran * 3 , but perspective will be gained if we discuss the more general case. First, observe in view of (12) that Second, whenever bounded, extracting a subsequence converging weakly to some ∈ 1,2 0 ( ) and using that the restriction map * 3 : 1,2 0 ( ) → 2 ( , R 3 ) is a fortiori continuous, by the Sobolev-Poincaré inequality (2), gives us in the limit that = * 3 ( ). To recap, we see from (12) that the quantity can be made arbitrarily small for appropriate , but if / ∈ * 3 [ 1,2 0 ( )] this is at the cost of letting ‖∇ ‖ 2 ( ,R 2 ) grow unbounded. In practice, one really measures So, given a priori bounds on ‖ ‖ 2 ( ,R 3 ) and ‖ ‖ 2 ( ) , 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): 3.1. Well posedness.
Proof. Since ‖∇ ‖ 2 ( ,R 2 ) is a norm equivalent to ‖ ‖ 1,2 ( ) on 1,2 0 ( ) by (3), the convex set is weakly compact in the Hilbert space 1,2 0 ( ). Then, since * 3 : is continuous hence also weakly continuous [8,Thm 3.10], the set of approximants 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 from the closed convex set , that can be written * 3 [ opt ] for some opt ∈ 1,2 0 ( ) which is unique because * 3 is injective. This ensures both existence and uniqueness of the solution opt to BEP.
Next, assume that ‖∇ opt ‖ 2 ( ,R 2 ) < . In this case, the minimum value of the criterion is achieved at opt which lies interior to the approximation set. We then and belong to ⊥ , see (12), thus . This contradicts the assumption on , thereby achieving the proof. (16) is not just closed but actually compact in 2 ( , 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.

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 onto the closed convex set 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 1,2 0 ( )-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.
i) A function opt ∈ 1,2 0 ( ) is the solution to BEP iff ‖∇ opt ‖ 2 ( ,R 2 ) = and there exists > 0 such that the following critical point equation holds, in the distributional sense on : ii) The function opt lies in ∞ ( ), and in the fractional Sobolev space 3/2,2 ( ), as well as in 1/2 ( ).
Proof. If we let : It is also instructive to establish (19) directly, using differentiation. For this, let ℳ ⊂ 1,2 0 ( ) be the smooth hypersurface comprised of those such that ‖∇ ‖ 2 ( ,R 2 ) = . The tangent space ℳ to ℳ at is the kernel of the linear form on 1,2 0 ( ) given by ↦ → ⟨∇ , ∇ ⟩ 2 ( ,R 2 ) . From Proposition 1 we know that opt ∈ ℳ and, by at opt must vanish on ℳ 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 ∈ 1,2 0 ( ), Clearly we get an equivalent equation upon restricting in (20) to range over ∞ 0 ( ), by density of the latter in 1,2 0 ( ). Thus, by the Green formula for Sobolev functions, we find that (20) is equivalent to ∀ ∈ ∞ 0 ( ), which means precisely that (19) holds in the distributional sense. To see that > 0, substitute = opt in (20) to obtain: and observe from (18), with = 0, that the above quantity is non-positive hence ≥ 0. Moreover ̸ = 0, otherwise (20) would imply that (17) is satisfied and then, arguing as in the proof of Proposition 1, it would entail that ∈ * Conversely, assume that ‖∇ opt ‖ 2 ( ,R 2 ) = and that (19) holds for some > 0. Then Equation (20) holds as well, and so does Equation (21) which is a special case of (20). Subtracting (21) from (20), we get that, ∀ ∈ 1,2 0 ( ), and by the Cauchy-Schwarz inequality we find that (18) holds when ‖∇ ‖ 2 ( ,R 2 ) ≤ . 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 1 2 2 . Since distributional derivatives commute, we see from (19) that ℎ := Δ opt a fortiori belongs to 2 ( ) because elements of the range of 3 are restrictions to of real analytic functions on R 2 . Now, if we let ( ) = − ∫︀ log | − |ℎ( ) d be the logarithmic potential of ℎ, we have that Δ = ℎ and it is a standard consequence of the Calderòn-Zygmund theory that ∈ 2,2 (R 2 ) [2, Thm 4.5.3]. Then, opt − is harmonic hence real analytic in , implying that the restriction of opt to any relatively compact disk in lies in 2,2 ( ). Since was arbitrary and 2,2 ( ) consists of continuous functions on by the Sobolev embedding theorem [12,Ch. 4], we get that opt ∈ ∞ ( ).

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 . For instance if is ∞ -smooth, then by elliptic regularity theory we get that opt is ∞ -smooth on , for Δ opt is real analytic there and opt vanishes on , see [21].
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 [3,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 = ‖∇ opt ‖ 2 ( ,R 2 ) when doing this. This is why we need to know more about the behavior of opt and as functions of . Hereafter, we use opt ( ) and ( ) = ‖∇ opt ( )‖ 2 ( ,R 2 ) when we want to make the dependence on explicit. 0 ( ) and R + respectively. As an illustration, let us mention that the curves obtained numerically in Section 5 and shown in Figures 2 and 3 provide a good overview of the properties listed in the lemma.
by the Cauchy-Schwarz inequality, the Green formula and the Poincaré-Sobolev inequality. Therefore, ‖Δ opt ‖ 2 ( ) cannot become small either. Because (19) implies that ) ⇒ ): assume that → 0 but nevertheless = ‖∇ opt ‖ 2 ( ,R 2 ) ̸ → ∞. Then, there exists a sequence → 0 for which ( ) is bounded, and we see from Equation (20) that converges weakly to in 2 ( , R 3 ): indeed, consider Γ ∈ 2 ( , R 3 ) and an arbitrary > 0 ; we decompose Γ as Γ = Γ ‖ + Γ ⊥ where where we applied that ‖ * 3 [ opt ] − ‖ 2 ( ,R 3 ) ≤ ‖ ‖ 2 ( ,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 * 3 [ opt ] − belongs to ⊥ (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 large enough (depending only on Γ and ), the last term is smaller than , which proves the announced weak convergence. Now, rearranging (21) as we find since 2 ( ) → 0 by assumption and ]︁ by hypothesis, therefore this limiting relation implies that ( ) → +∞, as indicated at the end of Section 2.3. This contradiction achieves the proof of the first assertion.
Next, we observe that for ∈ 1,2 0 ( ) with ̸ = 0, the Green formula and the fact that Δ = on 2 ( ) together imply: and setting = ( − 3  (23) and (24) that showing that 2 ( ) is smooth and strictly decreasing as a function of . Obviously the same holds for ( ), since it is never zero by the first assumption and the continuity of ↦ → ( ). Finally, since ( ) remains bounded as → +∞ by the first assertion of the lemma (or the monotonicity just proven), it follows from (2) that ‖ opt ‖ 2 ( ) remains bounded as well. Then (19) entails that Δ opt → 0 in 2 ( ) as → +∞. Applying the operator now shows that opt ( ) → 0 in 1,2 0 ( ) when → +∞, hence ( ) → 0. The monotonicity of ↦ → ( ) asserted in Lemma 4.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 (23) can be used for this purpose, however several issues arise in this connection. For instance the obvious boundary condition is opt ( ) → 0 in 1,2 0 ( ) 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 3 * 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 ′ an open set in R 2 such that ′ ⊂ , we define the local moment of the magnetization on ′ to be ⟨ ⟩ ′ = ∫︀ ′ . For magnetizations with a density like those considered in this paper (i.e. ∈ 2 ( , R 3 ) in our case), recovering the local moment on any disk ( , ) is tantamount to recovering because by Lebesgue's monotone convergence theorem ⟨ ⟩ ( , ) /( 2 ) → ( ) for a.e. ∈ when → 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 2 ( , R 3 ) with vanishing horizontal component are orthogonal to , by (12), we may set = (0, 0, ′ ) ∈ ⊥ in BEP and the results we obtained so far apply to any ′ . Thus, the numerical procedures described in sections to come can be used in principle to estimate the local moments of the vertical component of . This is consistent with the fact that compactly supported silent magnetizations have null vertical component.
The situation for the vector fields ( ′ , 0, 0) and (0, ′ , 0) is different, for they do not belong to ⊥ . However, if is an open neighborhood of ′ such that ⊂ , there exists a 1 -smooth function on which is identically 1 on ′ and 0 on ∖ . Then, the functions 1 ( 1 , 2 ) = ( 1 , 2 ) 1 and 2 ( 1 , 2 ) = ( 1 , 2 ) 2 are supported in and lie in 1,2 0 ( ), so the gradients ∇ 1 , ∇ 2 do lie in ⊥ and they coincide respectively with ( ′ , 0, 0) and (0, ′ , 0) on ′ . Now, if ′ does not intersect the support of , by compactness of the latter we can pick to contain no points of supp except those contained in ′ , and then ∫︀ ∇ . = ⟨ ⟩ ′ for = 1, 2. That is to say, if we know an open set ′ with ′ ⊂ whose boundary does not intersect supp (which means of course we have some a priori knowledge on ), we can construct an estimator of the local moment ⟨ ⟩ ′ .
In another connection, whereas ( ′ , 0, 0) and (0, ′ , 0) do not belong to ⊥ , their projections̃︀ 1, ′ and̃︀ 2, ′ onto that space in 2 ( , R 3 ) can be taken to be in BEP. Then, by (12) and the definition of the orthogonal projection, solving this BEP will provide us with an estimator of the local moment on ′ of the magnetization of minimum 2 ( , R 3 )-norm which produces the same field as . Note that̃︀ 1, ′ and︀ 2, ′ can be computed from the Hodge decomposition of ( ′ , 0, 0) and (0, ′ , 0) by solving a Neumann problem, identical to the one set up in the proof of [5,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. (19) has a unique solution opt ∈ 1,2 0 ( ) given ∈ 2 ( , R 3 ) and > 0. Indeed, taking the scalar product with an arbitrary ∈ 1,2 0 ( ), we may use the Green formula to obtain the equivalent equation:

Solving procedures. Equation
, where is the bound defined at Equation (10), we see that is continuous and coercive: where we can take, e.g., := /( 2 + 1) and := max( 2 , ), with the constant of the Sobolev-Poincaré inequality associated to (cf. Equation (2)). Therefore, appealing to the Lax-Milgram theorem [8,Cor. 5.8], Equation (25) admits a unique solution opt ∈ 1,2 0 ( ). Under the additional assumption that , Proposition 2 tells us that opt coincides with the solution of BEP corresponding to the level of constraint = ‖∇ opt ‖ 2 ( ,R 2 ) . A standard approach to numerically solve equations like (25) is to estimate the solutions in weak form, in a space of finite dimension. So, let ( ) ∈ be a finite family of functions 3 spanning a subspace = Span { , ∈ } of 1,2 0 ( ). Instead of Equation (25), we consider the restricted equation (26) ( where the solution is searched in the space : = ∑︀ ∈ for some real-valued coefficients ( ) ∈ . Considering Equation (26)   Solving this system yields the coefficients of the solution ∈ of Equation (26), which, by construction, is such that ( − opt , ) = 0 for any ∈ . Letting ‖ indicate the orthogonal projection of opt onto and ⊥ = opt − ‖ indicate its orthogonal complement, we thus get where the inequality is obtained by Cauchy-Schwarz theorem applied to the positive form . This reduces to | ( opt − , opt − )| 1/2 ≤ | ( ⊥ , ⊥ )| 1/2 , whence, with the coercivity of : Notice that all functions in satisfy ‖ ⊥ ‖ 1,2 0 ( ) ≤ ‖ − opt ‖ 1,2 0 ( ) by Pythagorean theorem. This irreducible error is inherent in the choice of the finite dimensional space for representing an approximate solution. Of course, one needs to ensure a priori that ‖ ⊥ ‖ 1,2 0 ( ) is small enough. Typically, the space will be controlled by some parameter such that the distance of any given function to = tends to zero when → 0. In our examples of Section 5, for instance, corresponds to the step size of the chosen mesh grid on .
In practice, the overall process is as follows. We fix a value > 0, so that the dimension of 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 , we iterate the following procedure, starting with some initial value = 0 > 0: compute the coefficients = ( ) for ∈ , together with the corresponding approximation to opt , and the associated constraint level ‖∇ ‖ 2 ( ,R 2 ) . If the latter is within a satisfactory range of , we stop. Otherwise we bisect with respect to variable , according to the rule indicated by Lemma 4.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 (i.e., ‖ ⊥ ‖ 1,2 0 ( ) 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 ‖ ⊥ ‖ 1,2 0 ( ) and the computational burden on dim( ). It is therefore important to pick a family ( ) 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 4.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 ( ) and did favor simplicity. More sophisticated choices of the basis are left here for further study.
The previous approach aims at approximating opt in 1,2 0 ( )-norm but tells us nothing about pointwise convergence, in particular it does not reflect that opt ∈ 1/2 ( ), see Proposition 2. This is why it seems worth describing another algorithm to solve Equation (19), 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 ∈ 2 ( , R 3 ), > 0, > 0, 1 ∈ 1,2 0 ( ) and ≥ 2, we consider the sequence of equations: 0 ( ) is the unknown and −1 acts as a parameter which was computed at the previous step.
Proof. Multiplying (19) by and subtracting it from (28), we obtain that , and taking the scalar product with − opt in 2 ( ) yields: Clearly we may assume that ̸ = opt for all ≥ 1, for if 0 = opt it follows from (30) that = opt for all ≥ 0 and there is nothing to prove. Now, using Equation (2) on the left-hand side and dividing on both sides by ‖ − opt ‖ we get Next, observe that ( being the bound introduced with Equation (10)) . Therefore, whenever 0 < < 1/ 2 , the operator − 3 * 3 is positive self-adjoint on 2 ( ), 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 ‖ − opt ‖ 2 ( ) decreases to 0 geometrically fast as → ∞.
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 has a smoothing effect, converting 2 ( )-convergence into 1/2 ( )-convergence. Moreover, if is a rectangle, the eigenfunctions and eigenvalues of the Laplacian are explicitly known (the eigenfunctions are products of sines) [18].
Therefore, once the products ⟨ * 3 [ , ], * 3 [ , ]⟩ 2 ( ,R 3 ) and ⟨ , , 3 [ ]⟩ 2 ( ) have been precomputed, it is easy to get the expansion of from the expansion of −1 . The computational burden in the case of a rectangle is overall fairly similar to the first approach, the precomputation of ⟨ * playing the same role as the precomputation of the matrix of Equation (27), and the precomputation of ⟨ , , 3 [ ]⟩ 2 ( ) corresponding to the computation of the right-hand side of Equation (26). However, thanks to the properties of , these expansions do not only converge in 2 ( ), but indeed in 1/2 ( ), while the first approach approximates opt only weakly. A thorough study of the precision required in the computation to make the 1/2 ( ) 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 1,2 0 ( ). Observe that Equation (19) is an elliptic partial differentialintegral equation to be solved on a square , subject to a homogeneous Dirichlet boundary condition. We thus make use of a family ( , ) of bilinear rectangular elements on a square mesh of . These elements are piecewise affine in the variables 1 , 2 separately, and are the simplest ones to expand 1,2 0 ( )-functions on a rectangle, see [10,Ch. 2] and [13].

Details of our implementation.
Mesh on . We choose a number ∈ N ⋆ that controls the number of elements that we use to mesh the rectangle . In order to simplify the presentation, we take to be the square [− , ] × [− , ] × {ℎ} ⊂ R 2 × {ℎ}, and we define = 2 /( + 1) and = − + for = 1, . . . , , so that the points of coordinates ( , ), 1 ≤ , ≤ constitute the × interior nodes of a square mesh of with step .
We denote with , the square centered at the node ( , ) with side-length equal to 2 and edges parallel to the axes. Observe that = ∪ , =1,··· , , but the union is not disjoint as the squares overlap (a generic , shares interior points with its 8 closest neighbors, ′ , ′ such that ′ ∈ { − 1, , + 1} and ′ ∈ { − 1, , + 1}). This is illustrated on Figure 1. Elements. We take for , = 1, · · · , the 1,2 0 ( )-functions , defined on by By construction, each , is 0 on the boundary and outside , (it actually belongs to Dot products on . We occasionally need to compute the dot product of 2 ( ) between a function of Span{ , } 1≤ , ≤ and a smooth function defined on . When it happens, we first approximate with the elements ( , ), i.e., on the component being considered), the value of ( ⋆̃︀)( , ) is evaluated, with the method described above, as the dot product . The numerical evaluations of , are obtained by an explicit formula for , analytically derived from its definition and Equation (6).
Right-hand side of the system. The right-hand side of Equation (26) requires to evaluate ⟨ , , 3 [ ]⟩ 2 ( ) for all 1 ≤ , ≤ . These are dot products on of the form described above, so we simply need to compute 3 [ ]( , ) for all 1 ≤ , ≤ . This is done using Equation (8): convolutions with the same kernels as before appear, except that the integrals involved are now on and not on . 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 = 0.00255 and ℎ = 0.00027, while the sample size is = 0.00197. This corresponds to a small, yet realistic, sample studied by SQUID microscopy. The discretization parameters and have been chosen to be = = 100.
We estimated the solution of the Critical Point Equation (19) using the abovedescribed algorithm, for different values of , and for the right-hand sides corresponding to = with = 1 and = 3. We denote with ( ) the approximation of the true solution opt so obtained. Notice that the case = 2 is similar to the case = 1, by symmetry (because 2 ( 1 , 2 ) = 1 ( 2 , 1 )), which is why we only show the results corresponding to 1 and 3 . Figures 2 and 3 show how the criterion ‖ * 3 [ ( )] − ‖ 2 ( ,R 3 ) and the level of constraint ( ) = ‖∇ ( )‖ 2 ( ,R 2 ) depend on . This is in accordance with the assertions of Lemma 4.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 [14]): when the level of constraint is fairly small, the criterion decreases very fast, meaning 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 intuitively means slightly larger oscillations of opt ). On the other hand, for larger values of , the criterion is barely improved even when 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 corresponding to = 10 −21 or = 10 −22 . We show on Figure 5 the functions ( ) ( = 1 and 3) for = 10 −21 . They are plotted on a rectangle slightly larger than to help understand how they behave at the boundary of (they are, by definition, equal to 0 outside ). 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 . Also, on the same bottom layer, the rectangles and are displayed, so as to recall their respective positions. We can see that the functions already have a fairly important oscillation pattern when = 10 −21 , more specifically on the region corresponding to ∖ . Indeed the maximal absolute value of 1 on ∖ is approximately 6.8 · 10 5 while its maximal absolute value on is roughly 0.89 · 10 5 .  (6)] ! e 3 k L 2 (S;R 3 ) ke 3 k L 2 (S;R 3 ) Figure 2. Approximation error of * 3 [ ( )] with respect to when varies (the scale for is logarithmic). The solid line corresponds to the case when = 1, while the dashed line corresponds to the case when = 3. As expected from Lemma 4.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 ( ) goes to 0, meaning that opt is forced to go to 0, whence the relative error tends to 1.
The same observation holds for 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 , and indeed the functions go to 0 fairly abruptly on . This means that, when estimating the moment ⟨ ⟩ from the data 3 [ ] 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 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{ , } 1≤ , ≤ and indicates that the choice of 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 ( ) with = 10 −21 to recover net moments, we considered a synthetic magnetization = ( 1 , 2 , 3 ) on . 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  (6) kr? e1 (6)k L 2 (Q;R 2 ) kr? e3 (6)k L 2 (Q;R 2 ) Figure 3. Constraint ( ) = ∇ ( ) as a function of (the scale for is logarithmic). The solid line corresponds to the case when = 1, while the dashed line corresponds to the case when = 3. As expected from Lemma 4.1, these are strictly decreasing smooth functions tending to +∞ when → 0 (i.e., log 10 ( ) → −∞) and tending to 0 when → +∞.
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 ⟨ ⟩ 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 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 3 .
We computed the values of 3 [ ] at the points of our mesh on . 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 3 [ ]. Using these values, we approximated the field and the noise as functions on using the elements ( , ), in order to evaluate the dot products with functions ( ) ( = 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 ⟨ ⟩ and their estimates decrease as tends to 0 because 3 [ ] is smooth enough that its integral against the linear estimators ( ) is fairly well evaluated. Interestingly, the component ⟨ 3 ⟩ is more accurately recovered than the others, which is likely connected with the fact that the inverse problem of recovering 3 from 3 [ ] has a unique solution (because silent sources have a null vertical component) whereas recovering 1 or 2 does not. This difference is reflected in the expression of the kernels generating 1 Λ and 2 Λ from , compare (5) and (7).  . The upper plot corresponds to the case when = 1, and the lower one to the case when = 3. As expected, the error decreases and tends to 0 as the constraint is relaxed.
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. 6. 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 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 * 3 [ , ] than those obtained using quadrature formulas, in order to get a more accurate version of the matrix and of the right-hand side in (26). 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 2 -density, and to more regular (e.g., Lipschitz-smooth) 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.