# Some Explicit Integrals Related to the Remanent Magnetization Inverse Problem

Sylvain Chevillard, Juliette Leblond

► To cite this version:

Sylvain Chevillard, Juliette Leblond. Some Explicit Integrals Related to the Remanent Magnetization Inverse Problem. RT-0525, Inria. 2025, pp.23. hal-05343567

HAL Id: hal-05343567

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

Submitted on 4 Nov 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
<figure><img src="image_2.png" /></figure>

# Some Explicit Integrals Related to the Remanent Magnetization Inverse Problem

Sylvain Chevillard, Juliette Leblond

## TECHNICAL REPORT

Nº 0525

November 2025

Project-Team Factas

ISSN 0249-0803 ISRN INRIA/RT-0525-FR+ENG
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

Some Explicit Integrals Related to the
Remanent Magnetization Inverse Problem

Sylvain Chevillard*, Juliette Leblond*

Project-Team Factas

Technical Report n° 0525 — November 2025 — 23 pages

**Abstract:** In this technical report, we recall explicit expressions of integrals related to the convolution of a polynomial with the Poisson kernel for which we also furnish simple proofs. Such integrals are ubiquitous in problems related to magnetic of gravity fields. In the specific case of the magnetization inverse problem for magnetic microscopy, we give a formula for the adjoint of the forward operator that maps the magnetization to the measured magnetic field. We give a completely explicit expression in the case when the adjoint operator is applied to a polynomial. In appendix, we provide a notebook of formulas for polynomials of small degree.

**Key-words:** explicit integrals, inverse magnetization problem, magnetic microscopy, QDM, SQUID, Poisson kernel, adjoint operator

<footnote>* Université Côte d'Azur, Inria, France</footnote>
Quelques intégrales explicites liées au problème inverse d'aimantation rémanente

**Résumé** : Dans ce rapport technique, nous rappelons les expressions explicites d'intégrales reliées à la convolution d'un polynôme et du noyau de Poisson pour lesquelles nous fournissons des preuves simples. Ces expressions sont omniprésentes dans les problèmes reliés aux champs magnétique ou de pesanteur. Pour le cas spécifique du problème inverse d'aimantation en microscopie magnétique, nous donnons une formule pour l'opérateur adjoint de l'opérateur direct qui à l'aimantation associe le champ magnétique mesuré. Nous donnons une expression complètement explicite dans le cas où l'opérateur adjoint est appliqué à un polynôme et fournissons en annexe un formulaire pour les polynômes de petits degrés.

**Mots-clés** : intégrales explicites, problème inverse d'aimantation, micro-
scopie magnétique, QDM, SQUID, noyau de Poisson, opérateur adjoint
# Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>2</td><td>An expression of the adjoint</td><td>5</td></tr><tr><td>3</td><td>Explicit form of the single integral</td><td>9</td></tr><tr><td>4</td><td>Explicit form of the double integral</td><td>11</td></tr><tr><td>4.1</td><td>Primitive of p(x)/(x<sup>2</sup> + a)<sup>1/2</sup></td><td>12</td></tr><tr><td>4.2</td><td>Primitive of p(x)/((x<sup>2</sup> + a)(x<sup>2</sup> + a + b)<sup>1/2</sup>)</td><td>13</td></tr><tr><td>5</td><td>References</td><td>15</td></tr><tr><td>A</td><td>Recall on inverse hyperbolic functions</td><td>16</td></tr><tr><td>B</td><td>Integrals expressed as sums of values of primitives</td><td>18</td></tr><tr><td>C</td><td>Value of the error bound</td><td>19</td></tr><tr><td>D</td><td>Collection of formulas</td><td>21</td></tr></table>

# 1 Introduction

The motivation of this report lies in the problem of recovering features of the remanent magnetization of rocks from measurements of the magnetic field that it generates on a horizontal plane above it. A typical use case that we target is magnetic microscopy where tiny samples of rocks are studied with an instrument such as a scanning superconducting quantum interference device (SQUID) [8] or a quantum diamond microscope (QDM) [5]. They provide a map of one or several components of the magnetic field on a horizontal rectangle just above the sample. Another possible application is the determination of magnetic anomalies in the ground from measurements of the field made from a plane overflying a region [9]. Our framework together with notations is recalled on Figure 1.

A rock sample assumed to be contained in a parallelepiped $A$ and carrying a magnetization density $m$ generates a magnetic field $B$ around it. Measurements are supposed to be performed on all points of a rectangle $Q$ lying at a height $x_3$ above the lower part of $A$. The rectangle $Q$ is not necessarily centered at the origin: its lower left corner is at coordinates $(\alpha_1, \alpha_2, x_3)$ while its upper right corner is at coordinates $(\beta_1, \beta_2, x_3)$. We assume that the measurement area does not intersect the magnetized sample. The particular (and more usual) case of a square with sides of half-size $R$ and centered at origin hence corresponds to $\beta_1 = \beta_2 = R = -\alpha_1 = -\alpha_2$.
<figure><img src="image_5.png" /><figcaption>Figure 1: Schematic setup with notations</figcaption></figure>

It has been long identified [6] that interesting information about the magnetization density can be retrieved from the values of integrals of the form $\iint_{\mathbb{R}^2} u(x_1, x_2) B_j(\mathbf{x}) dx_1 dx_2$ where $u$ is a smooth (and usually rather simple) function and $B_j(\mathbf{x})$ is the $j$-th component of the field, measured at point $\mathbf{x}$. In the case when $u$ is a polynomial of variables $x_1$ and $x_2$, explicit expressions of the integral over $Q$ (from which expressions of the integral over $\mathbb{R}^2$ can be obtained by taking the limit at infinity) are indeed known for long [7]. Though the existence of these exact expressions is acknowledged in the above cited work [6], the possibilities offered by accounting for the boundedness of the measurement area in the expression of the integral seem to have been neglected so far, while more or less sophisticated techniques were used to extrapolate the data outside $Q$ in order to properly evaluate the integral on the whole plane $\mathbb{R}^2$ [4].

The goal of the present report is to propose a comprehensive formulation of the adjoint of the magnetization-to-field operator, which is at the core of the integral methods described above. It develops in details computations that we have sketched in an extended abstract presented to a conference [3]. Specifically, we derive explicit formulas for the expression of the adjoint operator when applied to polynomials, with the hope that this report can be used as a convenient notebook for those who have to deal with this inverse magnetization problem. Though the explicit expressions of the involved integrals are easy to find in textbooks [10], it is not straightforward to find the path to derive them by integral calculus, and we thought it could be interesting for the reader to have a full description of it.
## 2 An expression of the adjoint

We start from the definition of the field: $B = -\mu_0 \nabla \phi$ where the scalar magnetic potential $\phi$ is given, for $x \notin A$, by

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

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

Let $j \in \{1, 2, 3\}$ denotes a component of the field. Accounting for the fact that $B_j$ is only known on $Q$, we identify $B_j$ with the forward operator which maps $\mathbf{m} \in [L^2(A)]^3$ to the corresponding component of the field $B_j[\mathbf{m}]$, seen as a function<fnref n="1" /> of $L^2(Q)$. Then, by definition, its adjoint operator $B_j^* : L^2(Q) \to [L^2(A)]^3$ satisfies, for all $u \in L^2(Q)$,

$$ \iint_Q u(x_1, x_2) B_j[\mathbf{m}](\mathbf{x}) \, dx_1 \, dx_2 = \iint_A \langle B_j^*[u](\mathbf{t}), \mathbf{m}(\mathbf{t}) \rangle \, d\mathbf{t}. \quad (2) $$

Let us denote with $S$ the left-hand side of Equation (2):

$$ S = \iint_Q u(x_1, x_2) B_j[\mathbf{m}](\mathbf{x}) \, dx_1 \, dx_2. $$

Starting from the definition of $S$, where $B_j[\mathbf{m}]$ is replaced by its definition and according to Equation (1), we get

$$ S = \frac{-\mu_0}{4\pi} \iint_Q u(x_1, x_2) \partial_{x_j} \iint_A \frac{\langle \mathbf{m}(\mathbf{t}), \mathbf{x} - \mathbf{t} \rangle}{\|\mathbf{x} - \mathbf{t}\|^3} d\mathbf{t} \, dx_1 \, dx_2. $$

Observing that $\partial_{x_j}$ commutes with $\iint_A$ (which is about integrating with respect to variable $\mathbf{t}$, not $x_j$) we obtain:

$$ S = \frac{-\mu_0}{4\pi} \iint_Q u(x_1, x_2) \iint_A \partial_{x_j} \left( \frac{\langle \mathbf{m}(\mathbf{t}), \mathbf{x} - \mathbf{t} \rangle}{\|\mathbf{x} - \mathbf{t}\|^3} \right) d\mathbf{t} \, dx_1 \, dx_2. $$

Now, since for any $(x_1, x_2)$, the quantity $u(x_1, x_2)$ is constant with respect to $\mathbf{t}$, we have:

$$ S = \frac{-\mu_0}{4\pi} \iint_Q \iint_A u(x_1, x_2) \partial_{x_j} \left( \frac{\langle \mathbf{m}(\mathbf{t}), \mathbf{x} - \mathbf{t} \rangle}{\|\mathbf{x} - \mathbf{t}\|^3} \right) d\mathbf{t} \, dx_1 \, dx_2. $$

Using Fubini's theorem we deduce:

$$ S = \iint_A \frac{-\mu_0}{4\pi} \iint_Q u(x_1, x_2) \partial_{x_j} \left( \frac{\langle \mathbf{m}(\mathbf{t}), \mathbf{x} - \mathbf{t} \rangle}{\|\mathbf{x} - \mathbf{t}\|^3} \right) dx_1 \, dx_2 \, d\mathbf{t}. $$

<footnote n="1">¹By abuse of notation, we identify here $Q = [\alpha_1, \beta_1] \times [\alpha_2, \beta_2] \times \{x_3\}$ (where $x_3$ is a fixed value) which is, really, a subset of $\mathbb{R}^3$, with $[\alpha_1, \beta_1] \times [\alpha_2, \beta_2] \subset \mathbb{R}^2$.</footnote>
Now, $\mathbf{m}(\mathbf{t})$ being constant with respect to $x_j$, it is a constant with respect to the $\partial_{x_j}$ operator, whence:

$$S = \iiint_{\mathcal{A}} \frac{-\mu_0}{4\pi} \iint_{Q} u(x_1, x_2) \left\langle \mathbf{m}(\mathbf{t}), \partial_{x_j} \left( \frac{\mathbf{x} - \mathbf{t}}{\|\mathbf{x} - \mathbf{t}\|^3} \right) \right\rangle \mathrm{d}x_1 \mathrm{d}x_2 \mathrm{d}\mathbf{t}.$$

By linearity of the inner product this can finally be rewritten as

$$S = \iiint_{\mathcal{A}} \left\langle \mathbf{m}(\mathbf{t}), \frac{-\mu_0}{4\pi} \iint_{Q} u(x_1, x_2) \partial_{x_j} \left( \frac{\mathbf{x} - \mathbf{t}}{\|\mathbf{x} - \mathbf{t}\|^3} \right) \right\rangle \mathrm{d}x_1 \mathrm{d}x_2 \mathrm{d}\mathbf{t}.$$

Let us denote with $\mathbf{B}_j^*[u]_{|i}$ the $i$-th component of $\mathbf{B}_j^*[u]$ (where $i \in \{1, 2, 3\}$). In view of Equation (2) and by identification, we deduce

$$\mathbf{B}_j^*[u]_{|i}(\mathbf{t}) = \frac{-\mu_0}{4\pi} \iint_Q u(x_1, x_2) \partial_{x_j} \left( \frac{x_i - t_i}{\|\mathbf{x} - \mathbf{t}\|^3} \right) \mathrm{d}x_1 \mathrm{d}x_2.$$

We do the change of variable $\mathbf{y} = \mathbf{x} - \mathbf{t}$ in the latter: for this purpose we introduce the notations $Q_t = \{(x_1 - t_1, x_2 - t_2), \text{ for } (x_1, x_2) \in Q\}$ and $s_t(y_1, y_2) = (t_1 + y_1, t_2 + y_2)$. The equation then rewrites:

$$\mathbf{B}_j^*[u]_{|i}(\mathbf{t}) = \frac{-\mu_0}{4\pi} \iint_{Q_t} (u \circ s_t)(y_1, y_2) \partial_{y_j} \left( \frac{y_i}{\|\mathbf{y}\|^3} \right) \mathrm{d}y_1 \mathrm{d}y_2. \quad (3)$$

A simple computation shows that for any $i \in \{1, 2, 3\}$

$$\partial_{y_i} \left( \frac{y_i}{\|\mathbf{y}\|^3} \right) = \frac{\|\mathbf{y}\|^2 - 3y_i^2}{\|\mathbf{y}\|^5}. \quad (4)$$

This equation and symmetry arguments obviously show that

$$\sum_{i=1}^{3} \partial_{y_i} \left( \frac{y_i}{\|\mathbf{y}\|^3} \right) = 0 \quad \text{and} \quad \partial_{y_j} \left( \frac{y_i}{\|\mathbf{y}\|^3} \right) = \partial_{y_i} \left( \frac{y_j}{\|\mathbf{y}\|^3} \right), \quad (5)$$

whence, together with Equation (3), we deduce that the matrix $(\mathbf{B}_j^*[u]_{|i})_{i,j}$ is symmetric with a vanishing trace. In particular, its third column can entirely be deduced from the two first ones. We hence further suppose without loss of generality that $j \in \{1, 2\}$. Let us define $k$ such that $\{j, k\} = \{1, 2\}$.

Notice that until now, we did not use any assumption about $Q$: the results obtained so far are valid whatever the shape of $Q$. Starting now, we use the fact that $Q = [\alpha_1, \beta_1] \times [\alpha_2, \beta_2]$, hence $Q_t = [\alpha_1 - t_1, \beta_1 - t_1] \times [\alpha_2 - t_2, \beta_2 - t_2]$. Integrating Equation (3) by parts yields

$$\frac{4\pi}{\mu_0} \mathbf{B}_j^*[u]_{|i}(\mathbf{t}) = - \int_{\alpha_k - t_k}^{\beta_k - t_k} \left[ \frac{y_i (u \circ s_t)(y_1, y_2)}{\|\mathbf{y}\|^3} \right]_{y_j = \alpha_j - t_j}^{\beta_j - t_j} \mathrm{d}y_k \\ + \iint_{Q_t} \frac{y_i (\partial_{x_j} (u) \circ s_t)(y_1, y_2)}{\|\mathbf{y}\|^3} \mathrm{d}y_1 \mathrm{d}y_2.$$
For a function $(y_1, y_2) \mapsto f(y_1, y_2)$, we denote with $I_j[f](y_1, y_2)$ a primitive of $f(y_1, y_2)/\|y\|^3$ with respect to variable $y_j$ ($j = 1, 2$), and we denote with $J[f](y_1, y_2)$ a double primitive of $f(y_1, y_2)/\|y\|^3$ with respect to variables $y_1$ and $y_2$. This leads us to our final result (coding the bounds of the intervals with $\varepsilon_1$ and $\varepsilon_2$ as detailed in Appendix B):

$$ 
\begin{aligned}
\mathbf{B}_j^*[u]_i(\mathbf{t}) &= \frac{-\mu_0}{4\pi} \left[ \left[ I_k[y_i(u \circ s_t)] - J[y_i(\partial_{x_j}(u) \circ s_t)] \right]_{\alpha_1-t_1}^{\beta_1-t_1} \right]_{\alpha_2-t_2}^{\beta_2-t_2} \quad (6) \\
&= \frac{-\mu_0}{4\pi} \sum_{\substack{\varepsilon_1 \in \{-1, 1\} \\ \varepsilon_2 \in \{-1, 1\}}} \varepsilon_1 \varepsilon_2 C_{i,j}[u]_{\substack{y_1 = \frac{\alpha_1+\beta_1}{2} + \varepsilon_1 \frac{\beta_1-\alpha_1}{2} - t_1, y_2 = \frac{\alpha_2+\beta_2}{2} + \varepsilon_2 \frac{\beta_2-\alpha_2}{2} - t_2 \\ \varepsilon_1 \in \{-1, 1\} \\ \varepsilon_2 \in \{-1, 1\}}}}
\end{aligned}
$$

where $C_{i,j}[u] = (I_k[y_i(u \circ s_t)] - J[y_i(\partial_{x_j}(u) \circ s_t)]) (y_1, y_2)$. The case when $u(x_1, x_2)$ is a polynomial is of particular interest: the problem boils down to integrating (once or twice) a polynomial over $\|y\|^3$. As we shall see, it turns out that such integrals admit explicit forms. But before going to them, let us specify three cases where this turns out useful:

1. Consider the case when $Q$ is a large symmetric square $Q = [-R, R] \times [-R, R]$ where $R \gg 1$. Then, we have

$$ \mathbf{B}_j^*[u]_i(\mathbf{t}) = -\frac{\mu_0}{4\pi} \sum_{\substack{\varepsilon_1 \in \{-1, 1\} \\ \varepsilon_2 \in \{-1, 1\}}} \varepsilon_1 \varepsilon_2 C_{i,j}[u]_{\substack{y_1 = \varepsilon_1 R - t_1, y_2 = \varepsilon_2 R - t_2, \\ \varepsilon_1 \in \{-1, 1\} \\ \varepsilon_2 \in \{-1, 1\}}} $$

from which an asymptotic expression as $R$ grows and of any order is easy to obtain. We previously derived such asymptotic expressions [2] with a less direct and more involved method.

2. Consider the case when $u$ is a finite element of $\mathbb{Q}_1$-type:

$$ u(x_1, x_2) = \left(1 - \frac{|x_1 - \kappa_1|}{\delta}\right) \left(1 - \frac{|x_2 - \kappa_2|}{\delta}\right), $$

for $(x_1, x_2) \in [\kappa_1 - \delta, \kappa_1 + \delta] \times [\kappa_2 - \delta, \kappa_2 + \delta]$, and 0 otherwise. Here, $\delta$ is the step size, and $(\kappa_1, \kappa_2)$ is the point of interpolation. We made use of such finite elements for the practical computation of a best linear evaluator in a previous work [1]: there, $\mathbf{B}_3^*[u](\mathbf{t})$ was approximated by a rough quadrature scheme. But we can actually derive an exact expression of it. Indeed, $u$ can be rewritten as the sum of four functions whose respective support are the corresponding quadrants of the elementary square, each of them being a polynomial of degree 2. Encoding the quadrant with the pair of variable $(\mu_1, \mu_2)$, where $\mu_1 \in \{-1, 1\}$ and $\mu_2 \in \{-1, 1\}$, we note

$$ u_{\mu_1, \mu_2}(x_1, x_2) = \left(1 - \mu_1 \frac{x_1 - \kappa_1}{\delta}\right) \left(1 - \mu_2 \frac{x_2 - \kappa_2}{\delta}\right), \quad (7) $$
<figure><img src="image_6.png" /><figcaption>Figure 2: Finite element of $\mathbb{Q}_1$-type</figcaption></figure>

defined on $[\alpha_1, \beta_1] \times [\alpha_2, \beta_2]$ where $[\alpha_n, \beta_n] = [\kappa_n - \delta, \kappa_n]$ if $\mu_n = -1$ and $[\alpha_n, \beta_n] = [\kappa_n, \kappa_n + \delta]$ if $\mu_n = 1$ (for $n \in \{1, 2\}$). We hence have $(\alpha_n + \beta_n)/2 = \kappa_n + \mu_n\delta/2$ and $(\beta_n - \alpha_n)/2 = \delta/2$, which leads us to the following explicit expression for $(-4\pi/\mu_0)\mathbf{B}_j^*[u]_i(t)$:

$$ \sum \varepsilon_1 \varepsilon_2 C_{i,j}[u_{\mu_1, \mu_2}]|_{y_1=\kappa_1+(\mu_1+\varepsilon_1)\frac{\delta}{2}-t_1, y_2=\kappa_2+(\mu_2+\varepsilon_2)\frac{\delta}{2}-t_2'} $$

where the sum is taken on the 16 possible values $(\varepsilon_1, \varepsilon_2, \mu_1, \mu_2) \in \{-1, 1\}^4$.

3. Consider the case when one needs to evaluate $\mathbf{B}_j^*[u]_i(t)$ for some sufficiently smooth (yet non-polynomial) $u$. One may consider approximating $u$ by a polynomial $p$. According to what we saw above, $\mathbf{B}_j^*[p]_i(t)$ admits the explicit expression of Equation (6) and the error

$$ |\mathbf{B}_j^*[p]_i(t) - \mathbf{B}_j^*[u]_i(t)| = |\mathbf{B}_j^*[p - u]_i(t)| $$

is easily bounded thanks to Equation (3):

$$ |\mathbf{B}_j^*[p-u]_i(t)| \le \frac{\mu_0}{4\pi} \|u-p\|_{\infty} \iint_{Q_t} \left| \partial_{y_j} \left( \frac{y_i}{\|\mathbf{y}\|^3} \right) \right| dy_1 dy_2, $$

where the supremum norm $\|u - p\|_{\infty}$ is taken on $Q$. The integral on $Q_t$ is upper-bounded by the integral on the whole $\mathbb{R}^2$. The function $(y_1, y_2, y_3) \mapsto y_i/\|\mathbf{y}\|^3$ is either even or odd with respect to any of its variable, and so are its partial derivatives. So the integral on the whole $\mathbb{R}^2$ is 4 times the integral on any semi-infinite quadrant. Therefore:

$$ |\mathbf{B}_j^*[p - u]_i(t)| \le \frac{\mu_0}{\pi} \|u - p\|_{\infty} \int_0^{+\infty} \int_0^{+\infty} \left| \partial_{y_j} \left( \frac{y_i}{\|\mathbf{y}\|^3} \right) \right| dy_1 dy_2. $$
This double integral can be computed explicitly thanks to the formulas presented in the next section. The computation (detailed in Appendix C) shows that

$$ \int_{0}^{+\infty} \int_{0}^{+\infty} \left| \partial_{y_j} \left( \frac{y_i}{\|\mathbf{y}\|^3} \right) \right| dy_1 dy_2 = \frac{\gamma}{y_3} \quad (8) $$

where $\gamma = 1$ if $i = j$ and $\gamma = 2\pi/(3\sqrt{3})$ otherwise.

## 3 Explicit form of the single integral

Let us start with the problem of a single integration: we show how to get an explicit expression of $I_k[p]$ for any polynomial $p$. To keep things simple, we forget the previous notations and we focus now on integrating, with respect to the real variable $x$, a polynomial $p(x)$ over $(x^2 + a)^{3/2}$ where $a > 0$.

**Case when $p(x) = x$**.

The case when $p(x) = x$ is easy as we recognize a derivative:

$$ \int \frac{x}{(x^2 + a)^{3/2}} dx = \int -\left(\frac{1}{2} 2x\right) (x^2 + a)^{-3/2} dx = \left[ \frac{-1}{(x^2 + a)^{1/2}} \right]. $$

**Case when $p(x) = x^2$**.

For the case when $p(x) = x^2$, we start with an integration by parts:

$$ \begin{aligned} \int \frac{x^2}{(x^2 + a)^{3/2}} dx &= \int x \frac{x}{(x^2 + a)^{3/2}} dx = \int x \left( \frac{-1}{(x^2 + a)^{1/2}} \right)' dx \\ &= \left[ \frac{-x}{(x^2 + a)^{1/2}} \right] + \int \frac{1}{(x^2 + a)^{1/2}} dx. \end{aligned} $$

For the integral of $(x^2 + a)^{-1/2}$, we recognize a known primitive:

$$ \begin{aligned} \int \frac{1}{(x^2 + a)^{1/2}} dx &= \int \frac{1}{\sqrt{a} \left( \left( \frac{x}{\sqrt{a}} \right)^2 + 1 \right)^{1/2}} dx \\ &= \frac{1}{\sqrt{a}} \int \frac{1}{(y^2 + 1)^{1/2}} \sqrt{a} dy \quad \text{where } y = \frac{x}{\sqrt{a}} \\ &= \left[ \arcsinh(y) \right] \\ &= \left[ \arcsinh \left( \frac{x}{\sqrt{a}} \right) \right]. \end{aligned} \quad (9) $$

Putting things together, we obtain

$$ \int \frac{x^2}{(x^2 + a)^{3/2}} dx = \left[ \frac{-x}{(x^2 + a)^{1/2}} + \arcsinh \left( \frac{x}{\sqrt{a}} \right) \right]. $$
Notice that, since inverse trigonometric functions can be expressed thanks to the ln function, other equivalent expressions exist for the above primitive (see Appendix A).

**General case.**

Indeed, the above formula for the case when $p(x) = x$ suggests us something: when deriving an expression of the form $x^m/(x^2 + a)^{1/2}$ we get a polynomial over $(x^2+a)^{3/2}$. More precisely, for any $m \ge 1$, we get the following expression for its derivative:

$$ \frac{mx^{m-1}}{(x^2 + a)^{1/2}} + \frac{-x^{m+1}}{(x^2 + a)^{3/2}} = \frac{mx^{m-1}(x^2 + a) - x^{m+1}}{(x^2 + a)^{3/2}} \\ = \frac{(m-1)x^{m+1} + amx^{m-1}}{(x^2 + a)^{3/2}}. \qquad (10) $$

In the particular case when $m = 1$, this reduces to $a/(x^2 + a)^{3/2}$, whence

$$ \int \frac{1}{(x^2 + a)^{3/2}} dx = \left[ \frac{x}{a(x^2 + a)^{1/2}} \right]. \qquad (11) $$

And, for $m \ge 2$, we deduce from Equation (10) the following recursion formula:

$$ \int \frac{x^{m+1}}{(x^2 + a)^{3/2}} dx = \frac{1}{m-1} \left[ \frac{x^m}{(x^2 + a)^{1/2}} \right] - \frac{ma}{m-1} \int \frac{x^{m-1}}{(x^2 + a)^{3/2}} dx. $$

Gathering the results obtained so far, this allows for the computation of a primitive in the case when $p(x)$ is any monomial. In principle, this is enough to compute a primitive for any polynomial $p(x)$ by linear combination thereof, but it turns out that the previous results can be summed up in a general recursion formula dealing with a generic polynomial $p(x)$ at once.

Specifically, if $p(x) = \sum_{k=0}^{m} \gamma_k x^k$ where $m \ge 2$, it holds:

$$ \int \frac{\sum_{k=0}^{m} \gamma_k x^k}{(x^2 + a)^{3/2}} dx = \left[ \gamma_0 \frac{x}{a(x^2 + a)^{1/2}} + \beta_{-1} \arcsin h \left( \frac{x}{\sqrt{a}} \right) + \frac{\sum_{k=0}^{m-1} \beta_k x^k}{(x^2 + a)^{1/2}} \right], $$

where $\beta_{m+1} = \beta_m = 0$, $\beta_k = \frac{\gamma_{k+1} - (k+2)\alpha\beta_{k+2}}{k-1}$ for $k \in [0, m-1] \setminus \{1\}$, $\beta_{-1} = \gamma_2 - 3\alpha\beta_3$ and $\beta_1 = -\beta_{-1}$. Even though it is required that $m \ge 2$ for the recursion formula to work, it is not required that $\gamma_m$ be non-zero, so the result can also be used for polynomials of degree less than 2, by letting $m = 2$ and $\gamma_2$ (and possibly $\gamma_1$) be null if necessary.

To prove the result, observe that (since, by definition, $\beta_m = \beta_{m+1} = 0$) the sum of the right-hand side remains the same when indexed up to $k = m+1$.
Then differentiating the right-hand side leads to

$$
\begin{align*}
& \frac{\gamma_0}{(x^2 + a)^{3/2}} + \frac{\beta_{-1}(x^2 + a)}{(x^2 + a)^{3/2}} + \frac{-\beta_0 x}{(x^2 + a)^{3/2}} + \sum_{k=1}^{m+1} \frac{(k-1)\beta_k x^{k+1} + ka\beta_k x^{k-1}}{(x^2 + a)^{3/2}} \\
&= \frac{1}{(x^2 + a)^{3/2}} \left( \gamma_0 + a\beta_{-1} - \beta_0 x + \beta_{-1} x^2 \right. \\
& \qquad \left. + \sum_{k=2}^{m+2} (k-2)\beta_{k-1} x^k + \sum_{k=0}^{m} (k+1)a\beta_{k+1} x^k \right) \\
&= \frac{1}{(x^2 + a)^{3/2}} \left( \gamma_0 + a\beta_{-1} - \beta_0 x + \beta_{-1} x^2 + a\beta_1 + 2a\beta_2 x \right. \\
& \qquad \left. + \sum_{k=2}^{m} ((k-2)\beta_{k-1} + (k+1)a\beta_{k+1}) x^k \right)
\end{align*}
$$

which simplifies to $\sum_{k=0}^{m} \frac{\gamma_k x^k}{(x^2+a)^{3/2}}$ thanks to the definitions of the $\beta_k$.

4 Explicit form of the double integral

Let $p(x_1, x_2)$ be a polynomial of two variables. In order to get an expression of $J[p]$, i.e., integrate twice (with respect to both variables) $p(x_1, x_2)/(x_1^2 + x_2^2 + x_3^2)^{3/2}$ where $x_3 \neq 0$ is some constant, we start by integrating, say, over $x_1$. The polynomial $p$ rewrites as a polynomial of variable $x_1$, whose coefficients $\gamma_k(x_2)$ are polynomials of variable $x_2$: $p(x_1, x_2) = \sum_{k=0}^{m} \gamma_k(x_2) x_1^k$ for some $m$ that we can assume greater or equal to 2 without loss of generality. Denoting $a(x_2) = x_2^2 + x_3^2$, and using the result obtained above, we get

$$
\int \frac{p(x_1, x_2)}{(x_1^2 + x_2^2 + x_3^2)^{3/2}} dx_1 = \frac{x_1 \gamma_0(x_2)}{a(x_2)(x_1^2 + x_2^2 + x_3^2)^{1/2}} \\
+ \beta_{-1}(x_2) \arcsin h \left( \frac{x_1}{\sqrt{a(x_2)}} \right) + \frac{\sum_{k=0}^{m-1} x_1^k \beta_k(x_2)}{(x_1^2 + x_2^2 + x_3^2)^{1/2}},
$$

where $\beta_{m+1} = \beta_m = 0$, $\beta_k(x_2) = \frac{\gamma_{k+1}(x_2)-(k+2)a(x_2)\beta_{k+2}(x_2)}{k-1}$ for $k \in [0, m-1] \setminus \{1\}$, and $\beta_{-1}(x_2) = \gamma_2(x_2) - 3a(x_2)\beta_3(x_2)$ and $\beta_1(x_2) = -\beta_{-1}(x_2)$. It is immediate from these recursion formulas that the $\beta_k(x_2)$ are indeed polynomials in variable $x_2$.

Now that it comes to integrate with respect to $x_2$, we actually have three
different expressions to integrate, namely:

$$ \frac{q_a(x_2)}{(x_2^2 + x_3^2)(x_2^2 + x_1^2 + x_3^2)^{1/2}}, \\ q_b(x_2) \operatorname{arcsinh}\left(\frac{x_1}{(x_2^2 + x_3^2)^{1/2}}\right), \\ \text{and } \frac{q_c(x_2)}{(x_2^2 + x_1^2 + x_3^2)^{1/2}}, $$

where $q_a$, $q_b$ and $q_c$ are polynomials in variable $x_2$ (the coefficients of $q_c$ being themselves polynomials in variable $x_1$).

The second case reduces to the first one thanks to an integration by parts. Indeed, deriving the arccsinh expression (that was first obtained by an integration with respect to variable $x_1$) with respect to variable $x_2$ leads a fraction similar to the first of the above three cases (this can be seen by explicitly doing the differentiation, but is actually the signature of an underlying structure; see Appendix A for details). More precisely, denoting with $I_2[q_b]$ a primitive of $q_b$, we have

$$ \int q_b(x_2) \operatorname{arcsinh} \left( \frac{x_1}{(x_2^2 + x_3^2)^{1/2}} \right) dx_2 = \left[ I_2[q_b](x_2) \operatorname{arcsinh} \left( \frac{x_1}{(x_2^2 + x_3^2)^{1/2}} \right) \right]_{x_2}^{x_1} + x_1 \int \frac{x_2 I_2[q_b](x_2)}{(x_2^2 + x_3^2)(x_2^2 + x_1^2 + x_3^2)^{1/2}} dx_2. $$

Therefore, we just have to study the first and third of the cases presented above. Again, to keep notations simple, we will further consider integration of functions of the real variable $x$, and we will denote with $p(x)$ a generic polynomial, and with $a$ and $b$ two real parameters such that $a > 0$ and $b \ge 0$.

## 4.1 Primitive of $p(x)/(x^2 + a)^{1/2}$

The case of a polynomial over $(x^2 + a)^{1/2}$ is done similarly to what we have done in Section 3.

**Case when $p(x) = 1$ and when $p(x) = x$**

The case when $p(x) = 1$ has indeed been already considered in Equation (9):

$$ \int \frac{1}{(x^2 + a)^{1/2}} dx = \left[ \operatorname{arcsinh} \left( \frac{x}{\sqrt{a}} \right) \right]. $$

For the case $p(x) = x$, we recognize a derivative:

$$ \int \frac{x}{(x^2 + a)^{1/2}} dx = \left[ (x^2 + a)^{1/2} \right]. $$
**General case.**

For $m \ge 1$, observe that the derivative of $x^m(x^2 + a)^{1/2}$ is

$$ 
\begin{aligned}
mx^{m-1}(x^2 + a)^{1/2} + \frac{x^{m+1}}{(x^2 + a)^{1/2}} &= \frac{mx^{m-1}(x^2 + a) + x^{m+1}}{(x^2 + a)^{1/2}} \\
&= \frac{amx^{m-1} + (m+1)x^{m+1}}{(x^2 + a)^{1/2}}.
\end{aligned}
$$

This provides us with the following recursion formula:

$$ \int \frac{x^{m+1}}{(x^2 + a)^{1/2}} dx = \frac{1}{m+1} \left[ x^m (x^2 + a)^{1/2} \right] - \frac{ma}{m+1} \int \frac{x^{m-1}}{(x^2 + a)^{1/2}} dx. $$

Again, it is more convenient from a computation point of view to use a recursion formula dealing with a generic polynomial $p(x)$ rather than separate monomials. In the present case, if $p(x) = \sum_{k=0}^{m} \gamma_k x^k$ where $m \ge 1$, it holds:

$$ \int \frac{\sum_{k=0}^{m} \gamma_k x^k}{(x^2 + a)^{1/2}} dx = \left[ \beta_{-1} \arcsinh \left( \frac{x}{\sqrt{a}} \right) + \left( \sum_{k=0}^{m-1} \beta_k x^k \right) (x^2 + a)^{1/2} \right], $$

where $\beta_{m+1} = \beta_m = 0$, $\beta_k = \frac{\gamma_{k+1} - (k+2)a\beta_{k+2}}{k+1}$ for $k \in [0, m-1]$, and $\beta_{-1} = \gamma_0 - a\beta_1$. Also, as before, the formula can be used with $m = 1$ and $\gamma_1 = 0$ to handle the case of polynomials of degree 0. The proof goes exactly as before.

## 4.2 Primitive of $p(x)/((x^2 + a)(x^2 + a + b)^{1/2})$

**Case when $p(x) = 1$**.

To integrate $1/((x^2 + a)(x^2 + a + b)^{1/2})$, we start by a change of variable to transform the term $(x^2+a+b)^{1/2}$ to $(\sinh(t)^2+1)^{1/2}$, namely $x = \sinh(t)\sqrt{a+b}$ which, after simplification, gives

$$ \int \frac{1}{(x^2 + a)(x^2 + a + b)^{1/2}} dx = \int \frac{1}{(a+b)\sinh(t)^2 + a} dt. $$

This calls for the change of variable $u = \tanh(t)$ for which $du = (1 - u^2) dt$ and $\sinh(t)^2 = u^2/(1 - u^2)$:

$$ \int \frac{1}{(a+b)\sinh(t)^2 + a} dt = \int \frac{1}{bu^2 + a} du. $$

Finally the change of variable $z = u\sqrt{b/a}$ leads us to a known derivative:

$$ \int \frac{1}{bu^2 + a} du = \frac{1}{\sqrt{ab}} \int \frac{dz}{1 + z^2} = \left[ \frac{\arctan(z)}{\sqrt{ab}} \right]. $$
Now, coming back to the original variable, we have on the one hand $z = u\sqrt{b/a} = \tanh(t)\sqrt{b/a}$ and on the other hand $x = \sinh(t)\sqrt{a+b}$, whence

$$z = \sqrt{\frac{b}{a}} \frac{\frac{x}{\sqrt{a+b}}}{\sqrt{1 + \frac{x^2}{a+b}}} = \sqrt{\frac{b}{a}} \frac{x}{\sqrt{x^2 + a + b}},$$

and therefore the result

$$\int \frac{1}{(x^2 + a)(x^2 + a + b)^{1/2}} dx = \left[ \frac{1}{\sqrt{ab}} \arctan \left( \frac{x\sqrt{b}}{\sqrt{a}\sqrt{x^2 + a + b}} \right) \right]. \quad (12)$$

**Case when $p(x) = x$**.

We start the same way when $p(x) = x$ with the change of variable $x = \sinh(t)\sqrt{a+b}$. This gives

$$\int \frac{x}{(x^2 + a)(x^2 + a + b)^{1/2}} dx = \int \frac{\sinh(t)\sqrt{a+b}}{(a+b)\sinh(t)^2 + a} dt.$$

Now, the natural change of variable in this case is $u = \cosh(t)$ for which $du = \sinh(t) dt$ and $\sinh(t)^2 = u^2 - 1$:

$$\int \frac{\sinh(t)\sqrt{a+b}}{(a+b)\sinh(t)^2 + a} dt = \int \frac{\sqrt{a+b}}{(a+b)u^2 - b} du$$

Finally, we do a change of variable to transform $(a+b)u^2 - b$ into $z^2 - 1$, namely $z = \frac{\sqrt{a+b}}{\sqrt{b}}u$:

$$\int \frac{\sqrt{a+b}}{(a+b)u^2 - b} du = \int \frac{dz}{\sqrt{b}(z^2 - 1)} = \frac{1}{2\sqrt{b}} \left[ \ln \left( \frac{z-1}{z+1} \right) \right].$$

Notice that we do not put absolute values inside $\ln$ because $z$ is $u$ multiplied by a number strictly larger than 1 and $u$ is larger or equal to 1, so $z > 1$. We come back to the original variable:

$$z = \frac{\sqrt{a+b}}{\sqrt{b}} \sqrt{1 + \frac{x^2}{a+b}} = \frac{\sqrt{x^2 + a + b}}{\sqrt{b}},$$

whence the result

$$\int \frac{x}{(x^2 + a)(x^2 + a + b)^{1/2}} dx = \left[ \frac{1}{2\sqrt{b}} \ln \left( \frac{\sqrt{x^2 + a + b} - \sqrt{b}}{\sqrt{x^2 + a + b} + \sqrt{b}} \right) \right]. \quad (13)$$

Another equivalent expression using arctanh exists for the above primitive (see Appendix A).
General case.

When $p(x)$ is a generic polynomial, its Euclidean division with $(x^2 + a)$ leads to $p(x) = q(x)(x^2 + a) + r(x)$ where $q(x)$ is a polynomial and $r(x)$ is a polynomial of degree less or equal to 1. Therefore

$$ \int \frac{p(x)}{(x^2 + a)(x^2 + a + b)^{1/2}} dx = \int \frac{q(x)}{(x^2 + a + b)^{1/2}} dx + \int \frac{r(x)}{(x^2 + a)(x^2 + a + b)^{1/2}} dx $$

which boils down to previously obtained results.

Notice that the polynomials $q$ and $r$ can easily be computed from $p$ using a simple recursion formula: namely, if $p(x) = \sum_{k=0}^{m} \gamma_k x^k$, then $r(x) = \mu_{-2} + \mu_{-1}x$ and $q(x) = \sum_{k=0}^{m-2} \mu_k x^k$ where the terms $(\mu_k)$ are defined thanks to $\mu_m = \mu_{m-1} = 0$ and $\mu_k = \gamma_{k+2} - a\mu_{k+2}$ for $k \in [-2, m-2]$. To see it, observe by a straightforward induction on the decreasing index $k$ that, for any $k \in [0, m-1]$, we have

$$ \left( \sum_{j=k}^{m-2} \mu_j x^j \right) (x^2 + a) = \left( \sum_{j=k+2}^{m} \gamma_j x^j \right) + a\mu_{k+1}x^{k+1} + a\mu_k x^k. $$

Now, defining $q(x) = \sum_{k=0}^{m-2} \mu_k x^k$ and $r(x) = \mu_{-2} + \mu_{-1}x$ and applying the above formula with $k=0$ we get

$$ q(x)(x^2 + a) = (p(x) - \gamma_1 x - \gamma_0) + a\mu_1 x + a\mu_0, $$

$$ \text{whence } p(x) - q(x)(x^2 + a) = (\gamma_1 - a\mu_1)x + (\gamma_0 - a\mu_0) = r(x). $$

5 References

[1] L. Baratchart, S. Chevillard, D. P. Hardin, J. Leblond, E. A. Lima, and J.-P. Marmorat. Magnetic moment estimation and bounded extremal problems. *Inverse Problems and Imaging*, 13(1), 2019. DOI: 10.3934/ipi.2019003.

[2] L. Baratchart, S. Chevillard, J. Leblond, E. A. Lima, and D. Ponomarev. Asymptotic method for estimating magnetic moments from field measurements on a planar grid. Working paper available on open archive, 2018. https://inria.hal.science/hal-01421157.

[3] S. Chevillard, J. Leblond, and A. Yousfi. On analytic formulas useful to solve the inverse magnetization problem. In *10ème Conférence Européenne sur les Méthodes Numériques en Electromagnétisme (Numelec 2024)*, Toulouse, France, 2024. https://inria.hal.science/hal-04905813.
[4] D. A. Clark. Methods for determining remanent and total magnetisations of magnetic sources – A review. *Exploration Geophysics*, 45(4):271–304, 2014. DOI: 10.1071/EG14013.

[5] D. R. Glenn, R. R. Fu, P. Kehayias, D. Le Sage, E. A. Lima, B. P. Weiss, and R. L. Walsworth. Micrometer-scale magnetic imaging of geological samples using a quantum diamond microscope. *Geochemistry, Geophysics, Geosystems*, 18(8):3254–3267, 2017. DOI: 10.1002/2017GC006946.

[6] 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. DOI: 10.23689/fidgeo-3081.

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

[8] J. R. Kirtley and J. P. Wikswo. Scanning SQUID microscopy. *Annual Review of Materials Science*, 29:117–148, 1999. DOI: 10.1146/annualrev.matsci.29.1.117.

[9] 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. DOI: 10.1186/BF03351848.

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

A Recall on inverse hyperbolic functions

Let $x$ and $y$ be real numbers and write $X = \exp(x)$. By definition $\arcsinh(y) = x$ is equivalent to $y = \sinh(x) = \frac{1}{2}(X - \frac{1}{X})$. It follows that $X^2 - 2yX - 1 = 0$. This equation of second order has discriminant $4y^2 + 4 > 0$ and hence admits two real solutions, one being positive and the other being negative (since their product equals $-1$). Since $X = \exp(x)$, it has to be the positive solution and hence $X = y + \sqrt{y^2 + 1}$. Consequently

$$ \arcsinh(y) = \ln \left( y + \sqrt{y^2 + 1} \right). $$

Now, since $(y + \sqrt{1 + y^2})(-y + \sqrt{1 + y^2}) = 1$, we have $\ln(y + \sqrt{1 + y^2}) = -\ln(-y + \sqrt{1 + y^2})$ whence a second form:

$$ \arcsinh(y) = -\ln \left( -y + \sqrt{y^2 + 1} \right). $$
Combining both expressions, we also get a nice symmetric form:

$$\arcsinh(y) = \frac{1}{2} \ln \left( \frac{y + \sqrt{1 + y^2}}{-y + \sqrt{1 + y^2}} \right).$$

Also, notice that $\arctanh(y) = x$ is equivalent to $y = \tanh(x) = \frac{X^2-1}{X^2+1}$, whence $X^2 = \frac{1+y}{1-y}$ and therefore $\arctanh(y) = \frac{1}{2}\ln(\frac{1+y}{1-y})$. Applied to $y = a/b$ this gives $\arctanh(a/b) = \frac{1}{2}\ln(\frac{b+a}{b-a})$ from which we get a fourth form:

$$\arcsinh(y) = \arctanh \left( \frac{y}{\sqrt{1 + y^2}} \right).$$

We previously established that $\int \frac{1}{(x^2 + a)^{1/2}} dx = \left[ \arcsinh \left( \frac{x}{\sqrt{a}} \right) \right]$. This can therefore also be written

$$\int \frac{1}{(x^2 + a)^{1/2}} dx = \left[ \frac{1}{2} \ln \left( \frac{x + \sqrt{x^2 + a}}{-x + \sqrt{x^2 + a}} \right) \right] = \left[ \operatorname{arctanh} \left( \frac{x}{(x^2 + a)^{1/2}} \right) \right].$$

The following expressions also hold, but beware that these primitives differ from the ones above by a constant:

$$\begin{gathered} \int \frac{1}{(x^2 + a)^{1/2}} dx = \left[ \ln \left( x + \sqrt{x^2 + a} \right) \right] \\ \text{and} \\ \int \frac{1}{(x^2 + a)^{1/2}} dx = \left[ -\ln \left( -x + \sqrt{x^2 + a} \right) \right]. \end{gathered}$$

Finally, let us observe how these results instantiate on our main problem:

$$\int \frac{1}{(y_1^2 + y_2^2 + y_3^2)^{1/2}} dy_1 = \operatorname{arctanh} \left( \frac{y_1}{(y_1^2 + y_2^2 + y_3^2)^{1/2}} \right)$$

and interestingly, applying Equation (13) and the rewriting rules described above:

$$\int \frac{y_1}{(y_1^2 + y_2^2)(y_1^2 + y_2^2 + y_3^2)^{1/2}} dy_1 = -\frac{1}{y_2} \operatorname{arctanh} \left( \frac{y_2}{(y_1^2 + y_2^2 + y_3^2)^{1/2}} \right).$$

This shows that the functions $L_j : (y_1, y_2) \mapsto \operatorname{arctanh} \left( \frac{y_j}{(y_1^2+y_2^2+y_3^2)^{1/2}} \right)$ (for $j \in \{1, 2\}$) are central in the explicit expressions of the primitives and double primitives of expressions of the form $p(y_1, y_2)/(y_1^2+y_2^2+y_3^2)^{3/2}$ where $p(y_1, y_2)$ is a polynomial. The following commutative diagram (that, alone, contains
Equations (9), (11), and (13)) sums it up:

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

where, as usual, $j$ and $k$ are such that $\{j, k\} = \{1, 2\}$.

Notice the following other interesting commutative diagram (containing Equations (11) and (12)):

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

where $A : (y_1, y_2) \mapsto \arctan\left(\frac{y_1 y_2}{y_3(y_1^2 + y_2^2 + y_3^2)^{1/2}}\right)$.

## B Integrals expressed as sums of values of primitives

If $x \mapsto F(x)$ is a function, the usual notation $[F(x)]_a^b$ means $F(b) - F(a)$. Coding the bounds of the segment $[a, b]$ with a parameter $\epsilon \in \{-1, 1\}$ ($-1$ coding for the left bound, and $1$ coding for the right bound), we observe that

$$[F]_a^b = \sum_{\epsilon \in \{-1, 1\}} \epsilon F\left(\frac{a+b}{2} + \epsilon \frac{b-a}{2}\right).$$

Let us now turn to the case of a function of two variables $(y_1, y_2) \mapsto f(y_1, y_2)$. We consider a double primitive $F$ of $f$ with respect to both variables, and we note $F_2 = \partial_{y_2} F$ (which is hence a primitive of $f$ with respect to $y_1$).
According to the above observation, we have

$$
\begin{align*}
& \int_{a_2}^{b_2} \int_{a_1}^{b_1} f(y_1, y_2) \, dy_1 \, dy_2 \\
&= \int_{a_2}^{b_2} \left[ F_2(y_1, y_2) \right]_{y_1=a_1}^{b_1} \, dy_2 \\
&= \int_{a_2}^{b_2} \sum_{\epsilon_1 \in \{-1, 1\}} \epsilon_1 F_2(y_1, y_2) \Big|_{y_1=\frac{a_1+b_1}{2}+\epsilon_1 \frac{b_1-a_1}{2}} \, dy_2 \\
&= \sum_{\epsilon_1 \in \{-1, 1\}} \epsilon_1 \left[ F(y_1, y_2) \Big|_{y_1=\frac{a_1+b_1}{2}+\epsilon_1 \frac{b_1-a_1}{2}} \right]_{y_2=a_2}^{b_2} \\
&= \sum_{\substack{\epsilon_1 \in \{-1, 1\} \\ \epsilon_2 \in \{-1, 1\}}} \epsilon_1 \epsilon_2 F(y_1, y_2) \Big|_{y_1=\frac{a_1+b_1}{2}+\epsilon_1 \frac{b_1-a_1}{2}, y_2=\frac{a_2+b_2}{2}+\epsilon_2 \frac{b_2-a_2}{2}}.
\end{align*}
$$

C Value of the error bound

We present here the detailed computation of Equation (8). Let us denote
with $E_{i,j}$ the quantity

$$
E_{i,j} = \int_{0}^{+\infty} \int_{0}^{+\infty} \left| \partial_{y_j} \left( \frac{y_i}{\|\mathbf{y}\|^3} \right) \right| dy_1 dy_2,
$$

for $i \in \{1, 2, 3\}$ and $j \in \{1, 2, 3\}$.

We first study the case $i = j = 3$: according to Equation (4) we have

$$
E_{3,3} = \int_{0}^{+\infty} \int_{0}^{+\infty} \left| \frac{y_1^2 + y_2^2 - 2y_3^2}{\|\mathbf{y}\|^5} \right| dy_1 dy_2 = \int_{0}^{\pi/2} \int_{0}^{+\infty} \left| \frac{r^2 - 2y_3^2}{(r^2 + y_3^2)^{5/2}} \right| r dr d\theta,
$$

where the last equality is obtained by a change to polar coordinates (y₁, y₂) =
(r cos(θ), r sin(θ)). The expression within the absolute value is negative for
r < y₃√₂ and positive for r > y₃√₂, whence

$$
E_{3,3} = \int_{0}^{\pi/2} \left( - \int_{0}^{y_3 \sqrt{2}} \frac{r(r^2 - 2y_3^2)}{(r^2 + y_3^2)^{5/2}} dr + \int_{y_3 \sqrt{2}}^{+\infty} \frac{r(r^2 - 2y_3^2)}{(r^2 + y_3^2)^{5/2}} dr \right) d\theta.
$$

Now, observe that $(r^2 - 2y_3^2) = (r^2 + y_3^2) - 3y_3^2$ whence

$$
\begin{align*}
\frac{r(r^2 - 2y_3^2)}{(r^2 + y_3^2)^{5/2}} &= \frac{r}{(r^2 + y_3^2)^{3/2}} - 3y_3^2 \frac{r}{(r^2 + y_3^2)^{5/2}} \\
&= \partial_r \left( \frac{-1}{(r^2 + y_3^2)^{1/2}} + y_3^2 \frac{1}{(r^2 + y_3^2)^{3/2}} \right) = \partial_r \left( \frac{-r^2}{(r^2 + y_3^2)^{3/2}} \right).
\end{align*}
$$
This finally gives $E_{3,3} = \int_{0}^{\pi/2} 2 \frac{2y_3^2}{(3y_3^2)^{3/2}} d\theta = \frac{2\pi}{3\sqrt{3}} \frac{1}{y_3}$.

Observing that $E_{i,j} = E_{j,i}$ by Equation (5), we can further assume that $j \in \{1, 2\}$. As we did before, we let $k \in \{1, 2\}$ be defined so that $\{j, k\} = \{1, 2\}$ and we now study the value of $E_{i,j}$ by studying separately the three cases $i = j, i = k$ and $i = 3$.

**The case when $i = j$:** from Equation (4), we have

$$ \partial_{y_j} \left( \frac{y_j}{\|\mathbf{y}\|^3} \right) = \frac{y_k^2 + y_3^2 - 2y_j^2}{\|\mathbf{y}\|^5}. $$

It is therefore positive for $y_j < \mu(y_k)$ and negative for $y_j > \mu(y_k)$, where $\mu(y_k) = \sqrt{\frac{y_k^2+y_3^2}{2}}$. Hence,

$$ 
\begin{aligned}
E_{j,j} &= \int_{0}^{+\infty} \left( \int_{0}^{\mu(y_k)} \partial_{y_j} \left( \frac{y_j}{\|\mathbf{y}\|^3} \right) dy_j - \int_{\mu(y_k)}^{+\infty} \partial_{y_j} \left( \frac{y_j}{\|\mathbf{y}\|^3} \right) dy_j \right) dy_k \\
&= \int_{0}^{+\infty} 2 \frac{\mu(y_k)}{(y_k^2 + y_3^2 + \mu(y_k)^2)^{3/2}} dy_k \\
&= \int_{0}^{+\infty} \frac{4}{3^{3/2}(y_k^2 + y_3^2)} dy_k \\
&= \frac{4}{3^{3/2}y_3^2} \int_{0}^{+\infty} \frac{1}{\left(\frac{y_k}{y_3}\right)^2 + 1} dy_k \\
&= \frac{4}{3\sqrt{3}y_3^2} \int_{0}^{+\infty} \frac{1}{s^2 + 1} y_3 ds \\
&= \frac{2\pi}{3\sqrt{3}} \frac{1}{y_3}.
\end{aligned}
$$

**The case when $i = k$:** the function $y_j \mapsto y_k/\|\mathbf{y}\|^3$ is obviously decreasing, so its derivative is negative. Therefore,

$$ 
\begin{aligned}
E_{k,j} &= - \int_{0}^{+\infty} \int_{0}^{+\infty} \partial_{y_j} \left( \frac{y_k}{\|\mathbf{y}\|^3} \right) dy_j dy_k \\
&= - \int_{0}^{+\infty} \left[ \frac{y_k}{\|\mathbf{y}\|^3} \right]_{y_k=0}^{+\infty} dy_k \\
&= \int_{0}^{+\infty} \frac{y_k}{(y_k^2 + y_3^2)^{3/2}} dy_k \\
&= \left[ \frac{-1}{(y_k^2 + y_3^2)^{1/2}} \right]_{y_k=0}^{+\infty} \\
&= \frac{1}{y_3}.
\end{aligned}
$$
The case when $i = 3$: similarly to the previous case, we have

$$ 
\begin{aligned}
E_{3,j} &= - \int_{0}^{+\infty} \left[ \frac{y_3}{\|\mathbf{y}\|^3} \right]_{y_k=0}^{+\infty} dy_k \\
&= \int_{0}^{+\infty} \frac{y_3}{(y_k^2 + y_3^2)^{3/2}} dy_k \\
&= \left[ \frac{y_3 y_k}{y_3^2 (y_k^2 + y_3^2)^{1/2}} \right]_{y_k=0}^{+\infty} \quad (\text{according to Equation (11)}) \\
&= \frac{1}{y_3}.
\end{aligned}
$$

## D Collection of formulas

We give here a collection of reference formulas that are easily obtained thanks to the results presented in Sections 3 and 4. To keep formulas short enough, we use the following notations, already introduced elsewhere in the document:

$$y_3 = x_3 - t_3,$$

$$a_j = y_j^2 + y_3^2 = y_j^2 + (x_3 - t_3)^2,$$

$$L_j = \operatorname{arctanh}\left(\frac{y_j}{\|\mathbf{y}\|}\right) = \operatorname{arctanh}\left(\frac{y_j}{(y_1^2 + y_2^2 + (x_3 - t_3)^2)^{1/2}}\right),$$

$$A = \arctan\left(\frac{y_1 y_2}{y_3 \|\mathbf{y}\|}\right) = \arctan\left(\frac{y_1 y_2}{(x_3 - t_3)(y_1^2 + y_2^2 + (x_3 - t_3)^2)^{1/2}}\right).$$

As before, let $i \in \{1, 2, 3\}$, let $j \in \{1, 2\}$ and let $k$ such that $\{k, j\} = \{1, 2\}$. Also $I_k[f]$ denotes a primitive of $f(y_1, y_2)/\|\mathbf{y}\|^3$ with respect to variable $y_k$ and $J[f]$ denotes a double primitive of $f(y_1, y_2)/\|\mathbf{y}\|^3$ with respect to variables $y_1$ and $y_2$.

### Simple primitives.

$$I_k[1] = \frac{y_k}{a_j \|\mathbf{y}\|},$$

$$I_k[y_3] = y_3 I_k[1] = \frac{y_k y_3}{a_j \|\mathbf{y}\|},$$

$$I_k[y_k] = \frac{-1}{\|\mathbf{y}\|},$$

$$I_k[y_j] = y_j I_k[1] = \frac{y_1 y_2}{a_j \|\mathbf{y}\|},$$

$$I_k[y_3 y_i] = y_3 I_k[y_i],$$

$$I_k[y_k^2] = L_k - \frac{y_k}{\|\mathbf{y}\|},$$
$$ 
\begin{aligned}
I_k[y_j^2] &= y_j^2 I_k[1] = \frac{y_k y_j^2}{a_j \|\mathbf{y}\|}, \\
I_k[y_1 y_2] &= y_j I_k[y_k] = \frac{-y_j}{\|\mathbf{y}\|}.
\end{aligned}
$$

Double primitives.

$$ 
\begin{aligned}
J[1] &= \frac{A}{y_3}, \\
J[y_3] &= y_3 J[1] = A, \\
J[y_j] &= -L_k, \\
J[y_3 y_i] &= y_3 J[y_i], \\
J[y_j^2] &= -y_3 A + y_k L_j, \\
J[y_k y_j] &= -\|\mathbf{y}\|.
\end{aligned}
$$

Expression of the adjoint.

In view of using Equation (6), we give below for various $u(x_1, x_2)$, $i$ and $j$ the explicit expression of $C_{i,j}[u] = (I_k[y_i(u \circ s_t)] - J[y_i(\partial_{x_j}(u) \circ s_t)])(y_1, y_2)$. Notice that (unsurprisingly since the matrix $(B_j^*[u]_i)_{i,j}$ is symmetric, see the comment after Equation (5)) one systematically observe that $C_{i,j}[u] = C_{j,i}[u]$.

$$ 
\begin{aligned}
C_{1,1}[1] &= \frac{y_1 y_2}{a_1 \|\mathbf{y}\|}, \\
C_{2,1}[1] &= \frac{-1}{\|\mathbf{y}\|}, \\
C_{3,1}[1] &= \frac{y_2 y_3}{a_1 \|\mathbf{y}\|}, \\
C_{1,2}[1] &= \frac{-1}{\|\mathbf{y}\|}, \\
C_{2,2}[1] &= \frac{y_1 y_2}{a_2 \|\mathbf{y}\|}, \\
C_{3,2}[1] &= \frac{y_1 y_3}{a_2 \|\mathbf{y}\|}.
\end{aligned}
$$

$$ 
\begin{aligned}
C_{1,1}[x_1] &= L_2 + \frac{y_1(y_1 + t_1)y_2}{a_1 \|\mathbf{y}\|}, & C_{1,1}[x_2] &= \frac{y_1}{\|\mathbf{y}\|} \left( \frac{t_2 y_2}{a_1} - 1 \right), \\
C_{2,1}[x_1] &= L_1 - \frac{t_1 + y_1}{\|\mathbf{y}\|}, & C_{2,1}[x_2] &= L_2 - \frac{t_2 + y_2}{\|\mathbf{y}\|}, \\
C_{3,1}[x_1] &= -A + \frac{y_2(y_1 + t_1)y_3}{a_1 \|\mathbf{y}\|}, & C_{3,1}[x_2] &= \frac{y_3}{\|\mathbf{y}\|} \left( \frac{t_2 y_2}{a_1} - 1 \right), \\
C_{1,2}[x_1] &= L_1 - \frac{t_1 + y_1}{\|\mathbf{y}\|}, & C_{1,2}[x_2] &= L_2 - \frac{t_2 + y_2}{\|\mathbf{y}\|},
\end{aligned}
$$
$$C_{2,2}[x_1] = \frac{y_2}{\|\mathbf{y}\|} \left( \frac{t_1 y_1}{a_2} - 1 \right), \quad C_{2,2}[x_2] = L_1 + \frac{y_1(y_2+t_2)y_2}{a_2 \|\mathbf{y}\|},$$

$$C_{3,2}[x_1] = \frac{y_3}{\|\mathbf{y}\|} \left( \frac{t_1 y_1}{a_2} - 1 \right), \quad C_{3,2}[x_2] = -A + \frac{y_1(y_2+t_2)y_3}{a_2 \|\mathbf{y}\|}.$$

$$C_{1,1}[x_1 x_2] = \|\mathbf{y}\| + t_2 L_2 + \frac{y_1(y_1+t_1)}{\|\mathbf{y}\|} \left( \frac{t_2 y_2}{a_1} - 1 \right),$$

$$C_{2,1}[x_1 x_2] = y_3 A + t_2 L_1 + t_1 L_2 - \frac{(y_1 + t_1)(y_2 + t_2)}{\|\mathbf{y}\|},$$

$$C_{3,1}[x_1 x_2] = -t_2 A + y_3 L_1 + \frac{y_3(y_1 + t_1)}{\|\mathbf{y}\|} \left( \frac{t_2 y_2}{a_1} - 1 \right),$$

$$C_{1,2}[x_1 x_2] = y_3 A + t_2 L_1 + t_1 L_2 - \frac{(y_1 + t_1)(y_2 + t_2)}{\|\mathbf{y}\|},$$

$$C_{2,2}[x_1 x_2] = \|\mathbf{y}\| + t_1 L_1 + \frac{y_2(y_2 + t_2)}{\|\mathbf{y}\|} \left( \frac{t_1 y_1}{a_2} - 1 \right),$$

$$C_{3,2}[x_1 x_2] = -t_1 A + y_3 L_2 + \frac{y_3(y_2 + t_2)}{\|\mathbf{y}\|} \left( \frac{t_1 y_1}{a_2} - 1 \right).$$

From those formulas, it is easy to compute the expression of $C_{i,j}[u_{\mu_1,\mu_2}]$ where $u_{\mu_1,\mu_2}$ was defined in Equation (7). However, the formulas start being too long to be explicitly written here.
Inria

RESEARCH CENTRE
Centre Inria d'Université Côte d'Azur

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-0803