# On analytic formulas useful to solve the inverse magnetization problem

Sylvain Chevillard, Juliette Leblond, Anass Yousfi

► To cite this version:

Sylvain Chevillard, Juliette Leblond, Anass Yousfi. On analytic formulas useful to solve the inverse magnetization problem. NUMELEC 2024 - 10ème Conférence Européenne sur les Méthodes Numériques en Electromagnétisme, Jul 2024, Toulouse, France. hal-04905813

HAL Id: hal-04905813

https://inria.hal.science/hal-04905813v1

Submitted on 22 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On analytic formulas useful to solve the inverse magnetization problem

Sylvain Chevillard*, Juliette Leblond* and Anass Yousfi†

*Factas team. Inria Research Center of Université Côte d'Azur, France

†Factas team. Université Côte d'Azur, France

Emails: firstname.lastname@inria.fr

**Abstract**—This work takes its roots in the inverse magnetization problem, where the magnetization properties of rocks are sought from measurements of the magnetic field they generate, taken on a horizontal plane at some height above the rocks. The Poisson kernel being at the core of the corresponding magnetic potential, double integrals of functions over the cube of a distance are ubiquitous in this context. We present explicit formulas for the particular case when the numerator of the integrand is a polynomial of two variables and show how these formulas can be useful to address the inverse magnetization problem.

## I. INTRODUCTION

The general problem motivating this work is the problem of recovering some characteristics of the magnetization of a rock specimen from measurements of the generated magnetic field taken on a horizontal plane above it. This setup is typical (though, of course, with very different scales), e.g., of measurements obtained from a flight over a region presenting a magnetic anomaly [8] or from measurements of the magnetic field of a small rock sample obtained from a scanning magnetic microscope [1]. We model this setup as in Fig. 1: a magnetic sample, enclosed in a bounded volume $A$, carries a magnetization density $m = (m_1, m_2, m_3)$ which produces a magnetic field $B$ measured at points $x = (x_1, x_2, x_3)$ of a square $Q$ at some fixed distance above. From Maxwell's equations in the static case [6, Sec. 5.9.C], $B = -\mu_0 \nabla \phi$ where the scalar magnetic potential $\phi$ is given by

$$ \phi(x) = \frac{1}{4\pi} \iiint_A \frac{\langle m(t), x - t \rangle}{\|x - t\|^3} dt, \quad (1) $$

where $\langle x, y \rangle = x_1 y_1 + x_2 y_2 + x_3 y_3$ denotes the usual scalar product and $\|x\| = (x_1^2 + x_2^2 + x_3^2)^{1/2}$.

<figure><img src="image_2.png" /><figcaption>Figure 1. Schematic setup with notations</figcaption></figure>

The full recovery of $m(t)$ for all $t \in A$ from the observed data is an ill-posed problem (it does not even have a unique solution [2]) but an interesting characteristic of the

magnetization, its net moment $\langle m \rangle = \iint_A m(t) dt$, is well defined from the data (though it remains unstable with respect to noise and calls for regularization techniques [1]). For this purpose, an old idea (see, e.g., [4]) is to integrate a component $B_j$ of $B$ (where $j \in \{1, 2, 3\}$) against simple functions $u(x_1, x_2)$. When $u(x_1, x_2)$ is a polynomial of small degree, and because he popularized them, these are called Helbig's integrals [5]. Accounting for the fact that $B_j$ is only known on $Q$, we identify $B_j$ with the forward operator which maps $m \in [L^2(A)]^3$ to the corresponding component of the field $B_j[m]$, seen as a function<fnref n="1" /> of $L^2(Q)$. Then, by definition, its adjoint operator $B_j^*[u] : L^2(Q) \to [L^2(A)]^3$ satisfies

$$ \iint_Q u(x_1, x_2) B_j[m](x) dx_1 dx_2 = \iint_A \langle B_j^*[u](t), m(t) \rangle dt. $$

Any piece of information on the structure of $B_j^*[u]$ could hence be useful to understand the effect of integrating $u$ against magnetic field data.

## II. RESULTS

First, we observe that the $i$-th component $B_j^*[u]_i$ of $B_j^*[u]$ (where $i \in \{1, 2, 3\}$) can be made more explicit. Replacing $B_j[m](x)$ in the above equation by $\partial_{x_j} \phi$ where $\phi$ is given by Eq. (1), then using Fubini's theorem to swap the integrals with respect to $x$ and $t$ and considering only the part in factor of $m_i$, we get a first expression for $B_j^*[u]_i$, from which, doing the change of variable $y = t - x$, we get

$$ B_j^*[u]_i(t) = \frac{-\mu_0}{4\pi} \iint_{Q_t} (u \circ s_t)(y_1, y_2) \partial_{y_j} \left[ \frac{y_i}{\|y\|^3} \right] dy_1 dy_2, $$
where $Q_t = [t_1 - R, t_1 + R] \times [t_2 - R, t_2 + R]$ and $s_t(y_1, y_2) = (t_1 - y_1, t_2 - y_2)$.

Now, observing that $\partial_{y_j}(y_i/\|y\|^3) = \partial_{y_i}(y_j/\|y\|^3)$ and $\sum_{k=1}^3 \partial_{y_k}(y_k/\|y\|^3) = 0$, we see that the matrix $(B_j^*[u]_i)_{i,j}$ is symmetric with a vanishing trace. We hence further suppose without loss of generality that $j \in \{1, 2\}$.

After integrating by part, one is reduced to integrate (w.r.t. a single variable) $y_i(u \circ s_t)(y_1, y_2)/\|y\|^3$ and integrate (w.r.t. both variables) $y_i[\partial_{y_j}(u \circ s_t)(y_1, y_2)]/\|y\|^3$. The case when $u(y_1, y_2)$ is a polynomial is of particular interest: the problem boils down to integrating (once or twice) a polynomial over $\|y\|^3$. For small degrees [7, Ch. 7] or specific cases [9, Ch. 1.2.43], explicit formulas are available in textbooks, but as far as we know, the following formulas,

<footnote n="1">¹By abuse of notation, we identify here $Q = [-R, R]^2 \times \{x_3\}$ (where $x_3$ is a fixed value) which is, really, a subset of $\mathbb{R}^3$, with $[-R, R]^2 \subset \mathbb{R}^2$.</footnote>
expressing the double integral with a simple recurrence on the coefficients of the polynomial $u$, do not exist in literature.

**Proposition 1.** Let $u(x_1, x_2) = \sum_{0 \le k+m \le p} \alpha_{k,m} x_1^k x_2^m$, where $\alpha_{k,m} \in \mathbb{R}$, be a polynomial of degree $p \in \mathbb{N}$. Then

$$ \iint \frac{u(x_1, x_2)}{\|\mathbf{x}\|^3} dx_1 dx_2 = a A(\mathbf{x}) \\ + \left(\sum_{k=0}^{p-1} u_k x_1^k\right) L_2(\mathbf{x}) + \left(\sum_{m=0}^{p-1} v_m x_2^m\right) L_1(\mathbf{x}) \\ + \left(\sum_{0 \le k+m \le p-2} w_{k,m} x_1^k x_2^m\right) \|\mathbf{x}\|, \quad (2) $$

where $L_i(\mathbf{x}) = \operatorname{argtanh}(x_i/\|\mathbf{x}\|)$ ($i \in \{1, 2\}$), $A(\mathbf{x}) = \arctan(x_1x_2/(x_3\|\mathbf{x}\|))$ and $\{(w_{k,m}), (u_k), (v_m), a\}$ are real numbers.

*Sketch of the proof.* Differentiating both sides of Eq. (2) with respect to $x_1$ and $x_2$, then multiplying by $\|\mathbf{x}\|^3$, we get<fnref n="2" />

$$ \begin{aligned} u(x_1, x_2) = & \sum_{2 \le k+m \le p} W_{k,m} x_1^k x_2^m \\ & + x_2^2 \sum_{k=1}^{p+2} (k+1) U_{k+1} x_1^k + x_1^2 \sum_{m=1}^{p+2} (m+1) V_{m+1} x_2^m \\ & + \sum_{k=1}^{p+2} U_k x_1^k + v_2 x_1^2 + \sum_{k=1}^{p+2} V_k x_2^k + u_2 x_2^2 \\ & + ax_3 + x_3^2 v_1 + x_3^2 u_1 + x_3^4 w_{1,1}, \end{aligned} \quad (3) $$

where

$$ \begin{aligned} W_{k,m} &= f_{k,m}(w_{k-1,m-1}, w_{k-3,m+1}, w_{k+1,m-3}, w_{k+1,m-1} \\ & \qquad w_{k-1,m+1}, w_{k+1,m+1}), \\ U_k &= g_k(u_{k+1}, u_{k-1}, w_{k-3,1}, w_{k-1,1}, w_{k+1,1}), \\ V_k &= h_k(v_{k+1}, v_{k-1}, w_{1,k-3}, w_{1,k-1}, w_{1,k+1}). \end{aligned} \quad (4) $$

such that $(f_{k,m})$, $(g_k)$ and $(h_k)$ are known $\mathbb{R}$-valued linear functions. Equating corresponding monomials in Equation (3) yields a square linear system of equations with $\{(w_{k,m}), (u_k), (v_m), a\}$ as the unknowns. Moreover, due to the specific structure of Eq. (4), an efficient recursive scheme can be implemented to solve this system instead of directly inverting the coefficient matrix. Namely, starting from $n = p-2$ and decreasing until $n = 1$, we are led at each step to the resolution of two tridiagonal systems to get the vector $z_n = (v_{n+1}, w_{0,n}, \dots, w_{k,n-k}, \dots, w_{n,0}, u_{n+1})$ (where $z_n = 0$ for $n > p-2$). On the other hand, $z_0 = (v_1, w_{0,0}, u_1)$, $z_{-1} = (v_0, u_0)$ and $z_{-2} = a$ are obtained via straightforward computations. To illustrate, if $n$ is even, we obtain a system for the odd-indexed unknowns:

$$ \begin{pmatrix} b_1 & c_1 & 0 & \cdots & 0 \\ a_3 & b_3 & c_3 & \cdots & 0 \\ 0 & a_5 & b_5 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & a_{n+3} & b_{n+3} \end{pmatrix} \begin{pmatrix} v_{n+1} \\ \vdots \\ w_{\text{odd,odd}} \\ \vdots \\ u_{n+1} \end{pmatrix} = \begin{pmatrix} d_1 \\ \vdots \\ d_{\text{odd}} \\ \vdots \\ d_{n+3} \end{pmatrix} $$

where the $(a_k)$, $(b_k)$ and $(c_k)$ are known integers, and the $(d_k)$'s are functions of $z_{n+2}$, $z_{n+4}$ and the $(\alpha_{k,m})$. A similar system permits retrieving the even-indexed parameters.

### III. ILLUSTRATION ON AN EXAMPLE

In order to illustrate a possible use of the proposed formulas, we consider an example inspired from [1]. Therein, a method is proposed to compute good linear estimators for the total net moment of a planar rock sample. The method involves evaluating $B_3^*[u](t)$, for $t \in A = [-s, s]^2 \times \{0\}$, where $u$ is an eigenfunction of the Laplacian on the square, which involves (among other computations) evaluating

$$ I(f_1, f_2) = \iint_{Q_t} f_1(x_1) f_2(x_2) \frac{1}{\|\mathbf{x}\|^3} dx_1 dx_2, $$

where $f_1(x_1) = \cos(\frac{n_1 \pi x_1}{2R})$ and $f_2(x_2) = \sin(\frac{n_2 \pi x_2}{2R})$ with $n_1$ and $n_2$ two integers, and $Q_t$ is defined as before as $Q_t = [t_1 - R, t_1 + R] \times [t_2 - R, t_2 + R]$.

For the sake of the example, say that $n_1 = 6$ and $n_2 = 4$, $R = 2.55e-3$, $s = 1.97e-3$ and $x_3 = 0.27e-3$. Using the dedicated Sollya tool [3] we can compute polynomials $p_i(x_i)$ ($i = 1, 2$) with binary64 coefficients such that $|p_i(x_i) - f_i(x_i)| \le 1e-5$ for all $x_i \in [-s - R, s + R]$. To achieve it, we choose $p_1$ even and of degree 28, while $p_2$ is odd and of degree 21. Replacing the $f_i$ by $p_i$ leads an error

$$ |I(f_1, f_2) - I(p_1, p_2)| \le \iint_{Q_t} \frac{\varepsilon}{\|\mathbf{x}\|^3} dx_1 dx_2 = \left[ \left[ \frac{\varepsilon A(\mathbf{x})}{x_3} \right] \right]_{Q_t} $$

where $\varepsilon = (2 + 1e-5) \times 1e-5$ and $A$ is the function defined in Proposition 1, such that $\partial_{x_1}\partial_{x_2}A(\mathbf{x}) = x_3/\|\mathbf{x}\|^3$. Defining $u(x_1, x_2) = p_1(x_1)p_2(x_2)$, the result of the previous section subsequently allows us to evaluate $I(p_1, p_2)$ at any point $t = (t_1, t_2, 0)$ by a simple evaluation of the indefinite double integral of $u(x_1, x_2)/\|\mathbf{x}\|^3$ at the bounds of $Q_t$.

### REFERENCES

[1] L. Baratchart, S. Chevillard, D. P. Hardin, et al. Magnetic moment estimation and bounded extremal problems. *Inverse Problems and Imaging*, 13(1), 2019.

[2] L. Baratchart, S. Chevillard, and J. Leblond. Silent and equivalent magnetic distributions on thin plates. In P. Jaming, A. Hartmann, K. Kellay, et al., editors, *Harmonic Analysis, Function Theory, Operator Theory, and their Applications*, Theta Series in Advanced Mathematics, pages 11–27. The Theta Foundation, 2017.

[3] S. Chevillard, M. Joldeş, and C. Lauter. Sollya: An Environment for the Development of Numerical Codes. In K. Fukuda, J. van der Hoeven, M. Joswig, and N. Takayama, editors, *Mathematical Software - ICMS 2010, volume 6327 of Lecture Notes in Computer Science*, pages 28–31, Heidelberg, Germany, 2010. Springer.

[4] D. A. Clark. Methods for determining remanent and total magnetisations of magnetic sources – A review. *Exploration Geophysics*, 45(4):271–304, 2014.

[5] K. Helbig. Some integrals of magnetic anomalies and their relation to the parameters of the disturbing body. *Zeitschrift fur Geophysik*, 29(2):83–96, 1963.

[6] J. D. Jackson. *Classical Electrodynamics*. John Wiley & Sons, 1999.

[7] K. Jung. *Schwerkraftverfahren in der angewandten Geophysik*. Akademische Verlagsgesellschaft, Geest & Portig K.-G., 1961.

[8] J. D. Phillips. Can we estimate total magnetization directions from aeromagnetic data using Helbig's integrals? *Earth, Planets and Space*, 57(8):681–689, 2005.

[9] A. Prudnikov, Y. A. Brychkov, and O. Marichev. *Integrals and Series*, volume 1. Taylor & Francis, 1986. Originally published in Russian in 1983 as *Integrale i ryady* by Nauka.

<footnote n="2">²We conventionally assume that parameters with indices exceeding the predefined bounds are equal to zero.</footnote>