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

HAL
open science

# A layer potential approach to inverse problems in brain imaging

Paul Asensio, Jean-Michel Badier, Juliette Leblond, Jean-Paul Marmorat, Masimba Nemaire

► To cite this version:

Paul Asensio, Jean-Michel Badier, Juliette Leblond, Jean-Paul Marmorat, Masimba Nemaire. A layer potential approach to inverse problems in brain imaging. Journal of Inverse and Ill-posed Problems, 2023, 10.1515/jiip-2023-0041. hal-03977724v2

HAL Id: hal-03977724

https://inria.hal.science/hal-03977724v2

Submitted on 20 Dec 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A layer potential approach to inverse problems in brain imaging

P Asensio¹, J-M Badier², J Leblond¹ and J-P Marmorat³, M Nemaire¹,⁴

December 20, 2023

## Abstract

We study the inverse source localisation problem using the electrical potential measured point-wise inside the head with stereo-ElectroEncephaloGraphy (sEEG), the electrical potential measured point-wise on the scalp with ElectroEncephaloGraphy (EEG) or the magnetic flux density measured point-wise outside the head with MagnetoEncephaloGraphy (MEG). We present a method that works on a wide range of models of primary currents, in particular we give details for primary currents that are assumed to be smooth vector-fields that are supported on and normally oriented to the grey/white matter interface. Irrespective of the data used we also solve the transmission problem of the electric potential associated with a recovered source hence we solve the cortical mapping problem. To ensure that the electric potential and normal currents are continuous in the head, the electric potential is expressed as a linear combination of double layer potentials and the magnetic flux density is expressed as a linear combination of single layer potentials. Numerically, we solve the problems on meshed surfaces of the grey/white matter interface, cortical surface, skull and scalp. A main feature of the numerical approach we take is that on the meshed surfaces we can compute the double and single layer potentials exactly and at arbitrary points. Because we explicitly study the transmission of the electric potential in head when using magnetic data, the coupling of electric and magnetic data in the source recovery problem is made explicit and shows the advantage of using simultaneous electric and magnetic data. We provide numerical examples of the source recovery and inverse cortical mapping using synthetic data.

*Keywords:* inverse problems, layer potentials, brain imaging

*AMS classification:* 31A25, 65N21, 65R32

## 1 Introduction

In this work, we look at the problem of identifying active regions of the brain using the electric potential and/or the magnetic flux density associated with brain activity.
This entails a lot into how the electric potential and magnetic flux density are trans-
mitted in the head. The problem of source localisation and the inverse transmission
problem have largely been solved separately. To solve the source localisation problem,
the transmission of the electric potential in the heads needs to be understood as the
observed electric potential and magnetic flux density are dependent on the secondary
currents due to the inhomogeneous electric conductivity of the head. So far the influ-
ence of the secondary currents has accounted for by either expressing the secondary
currents in terms of the source, this is the path taken in the software MNE, see [1], or
solving the inverse transmission problem first and then using the results in the source
localisation problem as is done in for example [2] and [3]. In either case, only the source
is explicitly recovered in source localisation. On the other hand, to solve the inverse
transmission problem, the fact that the electric potential is harmonic outside the sup-
port of the source is exploited and representation formulas of harmonic functions are
then used, this leads to formulations of the inverse transmission problem that do no
require knowledge on the source, see for example [4] and the references therein. In
this work we will take the view that these problems are different aspects of the same
problem and aim to solve them together. By solving the problems of source localisation
and inverse transmission, we provide context to the results we get and reinforce the
conclusions that are made from these results. For source localisation see for example
[2], and [1]. The problem of the transmission of the electric potential in the head has
been studied before notably in [4], [5] and [6] where the so-called symmetric boundary
elements method was employed. The symmetric boundary elements method uses the
single and double layer potentials, the normal derivative of the double layer potential
and the adjoint of the trace of the double layer potential, and an application of Galerkin
methods to solve the problem of the transmission of electromagnetic fields associated
with brain activity. Instead, we aim to use only the single and double layer potentials
with an application of boundary elements methods. In some cases, the application of
the boundary elements methods can be replaced with the method of fundamental solu-
tions, see for example [7], [8] for details on the theory and applications of the method
of fundamental solutions. We model brain activity, which may also be referred to as a
source or primary cerebral current, as a vector-field whose components are elements of
a Banach space supported on the grey/white matter interface, which is in accordance
with the neurophysiological structure of pyramidal neurons in the cerebral cortex, see
for example [9]. We give a special focus to $[W^{1/2,2}]^3$-vector-fields as sources as a way
to clearly demonstrate the use of layer potentials in this endeavour. We use the formu-
lations introduced by Geselowitz for the electric potential and magnetic flux density,
see for example [10], which results in expressions of the forward models of the electric
potential and the magnetic flux density that have strong relations with single and dou-
ble layer potentials. When considering simple approximations of the head geometry
such as a spherical head model, the forward models have explicit analytic expressions.
Analytic expressions are difficult to obtain when considering realistic head geometries
such as those that can be obtained via segmentation of MRI images, see for example
[11]. In this case appropriate discretisations of the brain structures and the vector-fields
are required. On any meshed surface we can obtain exact expressions for the associated
single and double layer potentials defined on the surface hence exact expressions for the
electric potential and magnetic flux density. This offers improved numerical accuracy and the versatility of applying the forward model of electric potential to both EEG and intra-cranial recordings as in sEEG. To further study how EEG, MEG and sEEG work we invite the reader to look into [12], [13] and [9], respectively. A point to note is that in MEG the magnetic flux density measured typically using Superconducting Quantum Interference Devices (SQUIDs) placed at a distance away from the head and around the head. Recently, Optically Pumped Magnetometers (OPMs) have been used to measure the magnetic flux density, see for example, [14].

When considering the inverse source localisation problem, we solve a Tikhonov regularised problem where we find the source that minimises a functional which involves the forward model. Hence, improved accuracy of the forward models in turn improves the source identification. In addition to source identification, the forward models that we employ can be co-opted to solve the so-called inverse cortical mapping problem, see [4], which is an inverse problem of the transmission of the electric potential within the head. In [4] a method based on the symmetric boundary elements method for solving the inverse cortical mapping problem was presented. The key difference between the method in [4] and the one we present here is that the unknown source is required in the computations. Further, the method we present is such that the inverse cortical mapping problem can be solved given either electric or magnetic data associated with brain activity. Coupling the problems of source localisation and cortical mapping should in principle improve the accuracy of the source localisation. The forward model for the magnetic flux density, see (13), provides a natural coupling of electric potential and magnetic flux density and we use this coupling to solve the inverse source localisation problem and inverse cortical mapping problem with simultaneous sEEG and MEG data which has interesting practical applications.

The paper is organised as follows. In Section 2, we give an overview on layer potentials and Maxwell's equation to motivate their later applications. In Section 3, we look at the forward models for the electric and magnetic potential in inhomogeneous domains, which forms the bedrock of this paper. In Section 4, we apply the forward models to $[W^{1/2,2}]^3$-vector-fields normally oriented to the grey/white matter interface to illustrate the use of the layer potentials. We will also look at how to build the discrete version of the problem which is largely applicable with minor changes when the source is taken from other Banach spaces. In Section 5, we discuss the inverse problems for sEEG, EEG and MEG together with the inverse cortical mapping problem. We show that a solution always exists and propose an algorithm to solve the problem. In Section 6, we present numerical examples of the algorithm with the sources taken as $[W^{1/2,2}]^3$-vector-fields and as dipoles to show the versatility of the algorithm and what we are aiming for. We provide a conclusion and outlook in Section 7.
# 2 Preliminaries

## 2.1 Function Spaces

We assume that $\Omega \subset \mathbb{R}^m$, $m \ge 3$, is a bounded open Lipschitz domain with connected boundary. Hence, $\mathbb{R}^m \setminus \partial\Omega$ has two connected components, namely $\Omega^-$ := int $\partial\Omega$ and $\Omega^+ := \text{ext } \partial\Omega$.

When $\Omega \subset \mathbb{R}^m$ is Lebesgue measurable and $1 \le p \le \infty$, we let $L^p(\Omega)$ be the space of (equivalence classes of a.e. coinciding) $\mathbb{R}$-valued measurable functions on $\Omega$ whose absolute value to the $p$-th power is integrable, with norm

$$ \|g\|_{L^p(\Omega)} = \left( \int_{\Omega} |g(y)|^p dy \right)^{1/p} \quad (\text{ess. sup}_{\Omega} |g| \text{ if } p = \infty). $$

Given a functional space $X$, we write $[X]^m$ for the corresponding space of vector-fields with $m$ components, each of which lies in $X$. For example, $[L^p(\Omega)]^m$ is the space of $\mathbb{R}^m$-valued vector fields $\mathbf{M}$ on $\Omega$ whose components belong to $L^p(\Omega)$, with norm

$$ \|\mathbf{M}\|_{[L^p(\Omega)]^m} = \left( \int_{\Omega} |\mathbf{M}|^p dy \right)^{\frac{1}{p}} \quad (\text{ess. sup}_{\Omega} |\mathbf{M}| \text{ if } p = \infty). $$

We let $W^{1,p}(\Omega)$ indicate the Sobolev space of functions lying in $L^p(\Omega)$ together with their first distributional derivatives. It is a Banach space with norm

$$ \|g\|_{W^{1,p}(\Omega)} = \left( \|g\|_{L^p(\Omega)}^p + \|\nabla g\|_{(L^p(\Omega))^n}^p \right)^{1/p}, $$

where $\nabla g = (\partial_1 g, \dots, \partial_n g)^t$ denotes the gradient of $g$ and $\partial_j g$ means the derivative with respect to the $j$-th variable.

Moreover, when $1 < p < \infty$, each $g \in W^{1,p}(\Omega)$ has a trace on $\partial\Omega$, say $\psi$, that lies in the fractional Sobolev space $W^{1-1/p,p}(\partial\Omega)$, where for $0 < s < 1$ we let

$$ \|\psi\|_{W^{s,p}(\partial\Omega)} := \|\psi\|_{L^p(\partial\Omega)} + \left( \int_{\partial\Omega} \int_{\partial\Omega} \frac{|\psi(x) - \psi(y)|^p}{|x - y|^{m-1+sp}} d\mathcal{H}(x)d\mathcal{H}(y) \right)^{1/p}, $$

where $\mathcal{H} = \mathcal{H}_{m-1}$ is the $(m-1)$-dimensional Haussdorf measure of $\partial\Omega$.

Given a Banach space $X$, we denote its dual by $X^*$. For $1 \le p < \infty$ and $1/p+1/q = 1$, the dual of $L^p(\Omega)$ is $L^q(\Omega)$, isometrically under the pairing $\langle f, g \rangle = \int_{\Omega} fg$, that is, $(L^p(\Omega))^* = L^q(\Omega)$. We let $W_0^{1,p}(\Omega)$ stand for the closure in $W^{1,p}(\Omega)$ of $C_c^\infty(\Omega)$, the space of infinitely differentiable functions with compact support in $\Omega$. For $1 < p < \infty$, the dual space of $W_0^{1,p}(\Omega)$ is denoted by $W^{-1,q}(\Omega)$, with $1/p + 1/q = 1$, that is, $(W_0^{1,p}(\Omega))^* = W^{-1,q}(\Omega)$. Similarly, for $1 < p < \infty$, we let $W^{-s,q}(\partial\Omega)$ indicate the fractional Sobolev space of negative order $-s$ and exponent $q$ which is the dual space of $W^{s,p}(\partial\Omega)$, $1/p + 1/q = 1$, that is, $(W^{s,p}(\partial\Omega))^* = W^{-s,q}(\partial\Omega)$.
## 2.2 Double and single layer potentials

The double and single layer potentials of a density $\psi$ defined on $\partial\Omega$, whose smoothness will be made precise later on, are defined by

$$K\psi(x) = \frac{1}{\omega_m} \int_{\partial\Omega} \psi(y) \frac{(x-y)}{|x-y|^m} \cdot \nu(y) d\mathcal{H}(y), \quad x \in \mathbb{R}^m \setminus \partial\Omega,$$

and

$$S\psi(x) = \frac{1}{(m-2)\omega_m} \int_{\partial\Omega} \psi(y) \frac{1}{|x-y|^{m-2}} d\mathcal{H}(y), \quad x \in \mathbb{R}^m \setminus \partial\Omega,$$

respectively, where $\omega_m$ indicates the area of the unit sphere in $\mathbb{R}^m$ and $\nu$ is the outward pointing unit normal to $\partial\Omega$. Note that $K\psi$ and $S\psi$ are harmonic in $\mathbb{R}^m \setminus \partial\Omega$. For $x \in \partial\Omega$, the double layer potential is defined as the singular integral

$$\begin{aligned} K\psi(x) &= p.v. \frac{1}{\omega_m} \int_{\partial\Omega} \psi(y) \frac{(x-y)}{|x-y|^m} \cdot \nu(y) d\mathcal{H}(y) \\ &= \lim_{\epsilon \to 0} \frac{1}{\omega_m} \int_{|x-y|>\epsilon} \psi(y) \frac{(x-y)}{|x-y|^m} \cdot \nu(y) d\mathcal{H}(y). \end{aligned}$$

As soon as $\psi \in L^p(\partial\Omega)$ for some $1 < p < \infty$, it follows from [15, Thm 1.10] that $K\psi$ exists a.e. in $\partial\Omega$ and belongs to $L^p(\partial\Omega)$. Furthermore, for almost every $y \in \partial\Omega$, $K\psi(x)$ converges to $(-\frac{1}{2}Id + K)\psi(y)$ (resp. to $(\frac{1}{2}Id + K)\psi(y)$), as $x \to y$ non-tangentially in $\Omega^-$ (resp. $\Omega^+$), see [15, Thm 1.10] for a precise statement and references. It follows that the non-tangential limits on $\partial\Omega$ of the double layer potential from inside and outside $\Omega$, differ by the density $\psi$ of the double layer potential. The normal derivatives of the double layer potential are continuous across $\partial\Omega$, see for example [16, Thm 6.13].

It is known that, $(-\frac{1}{2}Id + K) : W^{1/q,p}(\partial\Omega) \to W^{1/q,p}(\partial\Omega)$ is invertible for $p \in [\frac{3}{2}, 3]$ and $1/p + 1/q = 1$, see [17, Thm 8.1]. Hence, the double layer potential on $\Omega^-$ is a famous tool for solving the Dirichlet problem for the Laplace equation, which is to find $w : \Omega^- \to \mathbb{R}$ such that

$$\begin{aligned} \Delta w &= 0 \text{ in } \Omega^-, \\ w &= g \text{ on } \partial\Omega, \end{aligned} \quad (1)$$

with $g \in W^{1/q,p}(\partial\Omega)$, $p \in [\frac{3}{2}, 3]$, for which the solution is $w = K(-\frac{1}{2}Id + K)^{-1}g$ and belongs to $W^{1,p}(\Omega^-)$. Here, the boundary condition in (1) is satisfied both as a Sobolev trace and as a non-tangential limit a.e.

Further, it is also known that $(\frac{1}{2}Id + K) : W^{1/q,p}(\partial\Omega)/\langle 1 \rangle \to W^{1/q,p}(\partial\Omega)/\langle 1 \rangle$ is invertible for $p \in [\frac{3}{2}, 3]$, where the quotient by $\langle 1 \rangle$ means “modulo constants”, see [17, Thm 8.1]. Likewise, the double layer potential on $\Omega^+$ becomes a tool for solving the exterior Dirichlet problem, which is to find $w : \Omega^+ \to \mathbb{R}$ such that

$$\begin{aligned} \Delta w &= 0 \text{ in } \mathbb{R}^m \setminus \overline{\Omega^-} \cup \{\infty\}, \\ w &= g \text{ on } \partial\Omega, \end{aligned}$$

with $g \in W^{1/q,p}(\partial\Omega)/\langle 1 \rangle$, $p \in [\frac{3}{2}, 3]$ thus $w = K(\frac{1}{2}Id + K)^{-1}g$ will solve the exterior Dirichlet problem up to a constant. Further, $w|_{\Omega^+ \cap B(0,R)} \in W^{1,p}(\Omega^+ \cap B(0, R))$ for all
$R > 0$ and $\nabla w \in [L^p(\Omega^+)]^m$. Here, as $m \ge 3$, harmonicity at infinity means that $\lim_{|x| \to \infty} w(x) = 0$ [18, Thm 4.8].

Unlike the double layer potential, the single layer potential is continuous across $\partial\Omega$, see for example [16, Thm 6.13], though its normal derivative is not. Note that the gradient of the single layer potential is given by

$$ \nabla S\psi(x) = -\frac{1}{\omega_m} \int_{\partial\Omega} \frac{(x-y)}{|x-y|^m} \psi(y) d\mathcal{H}(y), \quad x \in \mathbb{R}^m \setminus \partial\Omega. $$

For $1 < p < \infty$ and $\psi \in L^p(\partial\Omega)$, it holds for a.e. $y \in \partial\Omega$ that $\nu(y) \cdot \nabla S\psi(x)$ converges to $(\frac{1}{2}Id - K^*)\psi(y)$ (resp. $(-\frac{1}{2}Id + K^*)\psi(y)$) as $x \to y$ non-tangentially in $\Omega^-$ (resp. $\Omega^+$), where $K^*$ is the adjoint of $K$ as operators on $L^p(\partial\Omega)$:

$$ K^*\psi(x) = p.v. - \frac{1}{\omega_m} \int_{\partial\Omega} \psi(y) \frac{(x-y)}{|x-y|^m} \cdot \nu(x) d\mathcal{H}(y) $$

see [15, Thm 1.11] for a statement and further references.

It follows from [17, Thm 8.1] that $(\pm \frac{1}{2}Id + K^*) : \tilde{W}^{-1/p,p}(\partial\Omega) \to \tilde{W}^{-1/p,p}(\partial\Omega)$ extends to an invertible map for $p \in [\frac{3}{2}, 3]$, where we have set

$$ \tilde{W}^{-\frac{1}{p},p}(\partial\Omega) := \{ f \in W^{-\frac{1}{p},p}(\partial\Omega) : \int_{\partial\Omega} f(y) d\mathcal{H}(y) = 0 \}. $$

For the appropriate range of $p$, the single layer potential on $\Omega^-$ solves the Neumann problem for the Laplace equation

$$ \begin{gathered} \Delta w = 0 \text{ in } \Omega^{-} \\ , \nabla w \cdot \nu = g \text{ on } \partial\Omega, \end{gathered} \qquad (2) $$

where $g \in \tilde{W}^{-1/p,p}(\partial\Omega)$. Likewise, the single layer potential on $\Omega^+$ can be used to solve the exterior Neumann problem:

$$ \begin{gathered} \Delta w = 0 \text{ in } \mathbb{R}^m \setminus \overline{\Omega^{-}} \cup \{\infty\}, \\ \nabla w \cdot \nu = g \text{ on } \partial\Omega, \end{gathered} \qquad (3) $$

where $g \in \tilde{W}^{-1/p,p}(\partial\Omega)$. Thus, by [17, Thm 9.2], the solution to (2 (resp. (3) can be written as $w = S(\pm \frac{1}{2}Id + K^*)^{-1}g$, up to an additive constant. Moreover, $w$ belongs to $W^{1,p}(\Omega^-)$ (resp. $w|_{\Omega^+ \cap B(0,R)} \in W^{1,p}(\Omega^+ \cap B(0, R))$ for all $R > 0$ and $\nabla w \in [L^p(\Omega^+)]^m$).

## 2.3 Maxwell's equations

Henceforth, we take $m = 3$. According to [19], Maxwell's equations are given in differential form as follows

$$ \nabla \times \mathbf{E} = -\mu \frac{\partial \mathbf{B}}{\partial t}, \quad \nabla \times \mathbf{B} = \mathbf{J} + \epsilon \frac{\partial \mathbf{E}}{\partial t}, \quad \epsilon \nabla \cdot \mathbf{E} = \rho, \quad \mu \nabla \cdot \mathbf{B} = 0, $$

where $\mathbf{E} \in \mathbb{R}^3$ is the electric field, $\mathbf{B} \in \mathbb{R}^3$ is the magnetic flux density both of which are produced by the charge and current densities, $\rho \in \mathbb{R}$ and $\mathbf{J} \in \mathbb{R}^3$, respectively, $\mu \in \mathbb{R}_+$ and $\epsilon \in \mathbb{R}_+$ are the magnetic permeability and electric permittivity, respectively.
If the partial derivatives with respect to time in the above expressions are negligible then we can use a quasi-static approximation of the Maxwell equations given as follows

$$ \nabla \times \mathbf{E} = 0, \quad \nabla \times \mathbf{B} = \mathbf{J}, $$

For EEG and MEG the quasi-static approximation can be made for the electromagnetic dynamics. In the quasi-static regime $\nabla \times \mathbf{E} = 0$, thus it follows that $\mathbf{E} = -\nabla\phi$ for some scalar potential $\phi$. The total current density is of the form $\mathbf{J} = \mathbf{J}^i - \sigma\nabla\phi$ where $\mathbf{J}^i$ is the impressed primary current, $-\sigma\nabla\phi$ is the ohmic current and $\sigma$ is the electric conductivity which needs not be uniform in the domain. The elliptic equation

$$ \nabla \cdot (\sigma \nabla \phi) = \nabla \cdot \mathbf{J}^i, \qquad (4) $$

can be derived from the quasi-static approximation and we shall use it in this work.

# 3 Forward Models

## 3.1 Unbounded homogeneous domain

In much of the literature on inverse problems of source localisation the “elementary electromagnetic object” that is sort after is a current dipole. However, in this work we will suppose that the source is an element of $[\Xi]^3$ for some Banach space $\Xi$. Thus, it follows from (4) that the electrical potential associated with primary currents $\xi \in [\Xi]^3$ satisfies

$$ \nabla \cdot (\sigma \nabla \phi) = \nabla \cdot \boldsymbol{\xi}. \qquad (5) $$

From (5) in an infinite homogeneous medium the potential, $\phi$, is given by

$$ \sigma\phi(\boldsymbol{\xi})(x) = \frac{1}{4\pi} \int \frac{\nabla \cdot \boldsymbol{\xi}(y)}{|x - y|} dy, \quad x \notin \text{supp} \boldsymbol{\xi}, \qquad (6) $$

and the gradient of the potential, $\nabla\phi$, is given by

$$ \sigma \nabla \phi(\boldsymbol{\xi})(x) = \frac{1}{4\pi} \int \frac{\boldsymbol{\xi}(y)}{|x - y|^3} dy - 3 \int (x - y) \frac{(x - y)}{|x - y|^5} \cdot \boldsymbol{\xi}(y) dy, \quad x \notin \text{supp} \boldsymbol{\xi}. $$

Since the magnetic flux density, $\mathbf{B}$, is divergence free it is the curl of the vector magnetic potential which we denote by $\mathbf{A}$, that is,

$$ \mathbf{B} = -\nabla \times \mathbf{A}. $$

Given a primary current $\boldsymbol{\xi}$ then

$$ \mathbf{A}(\boldsymbol{\xi})(x) = \frac{\mu}{4\pi} \int \frac{\boldsymbol{\xi}(y)}{|x - y|} dy, \quad x \notin \text{supp} \boldsymbol{\xi}. \qquad (7) $$

Using the relationship above between the vector magnetic potential and the magnetic flux density we have that

$$ \mathbf{B}(\boldsymbol{\xi})(x) = \frac{\mu}{4\pi} \int \boldsymbol{\xi}(y) \times \nabla_y \frac{1}{|x - y|} dy, \quad x \notin \text{supp} \boldsymbol{\xi}. \qquad (8) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: An example cross-section of a non-homogeneous domain and the placement of intra-cranial electrodes for sEEG, in red; the blue arrows represent the normal orientation of the source term on its support.</figcaption></figure>

## 3.2 Non-homogeneous bounded domains

We now study the transmission of the electric potential and magnetic flux density in a bounded non-homogeneous conductor. Figure 1 is an example of such a conductor. To that end we consider the following bounded domain; let $\Omega \subset \mathbb{R}^3$ be a nested non-homogeneous bounded Lipschitz domain such that $\Omega_0, \Omega_1, \dots, \Omega_n$ are nested Lipschitz domains with $\Omega_0 \subset \Omega_1 \subset \dots \subset \Omega_n$ and $\Omega = \cup_i \Omega_i$. We let the boundaries $\partial\Omega_{i-1} \cap \partial\Omega_i = \Sigma_i$ with $\nu_i, i = 1, 2, \dots, n$ being the outward pointing unit normal to $\Sigma_i, i = 1, 2, \dots, n$, respectively. For consistency of notation the outer boundary of $\Omega_n$ shall be called $\Sigma_{n+1}$ with $\nu_{n+1}$ as the outward pointing unit normal to $\Sigma_{n+1}$. We will assume that the support of $\xi \in [\Xi]^3$ is a proper subset of $\bar{\Omega}_0$. The electric conductivities of the different domains are constant in each domain but different between domains, we call them $\sigma_0, \sigma_1, \dots, \sigma_n$ for $\Omega_0, \Omega_1, \dots, \Omega_n$, respectively. The conductivity outside $\Sigma_{n+1}$ will be set to zero, hence $\sigma_{n+1} = 0$. On any of the interfaces $\Sigma_i, i = 1, 2, \dots, n+1$, let $\sigma_i^-$ and $\sigma_i^+$ be the conductivities inside and outside, respectively.

The electric potential at the interfaces is denoted by $\phi_i^-$ and $\phi_i^+$ depending on whether the electric potential is taken as a non-tangential limit approaching the interface $\Sigma_i$ from inside or outside, respectively. On each interface $\Sigma_i$ the electric potential satisfies

$$ \phi_i^- = \phi_i^+, \\ \sigma_i^- \partial_{\nu_i} \phi^- = \sigma_i^+ \partial_{\nu_i} \phi^+, \qquad (9) $$

since $\phi$ is a solution to the elliptic problem (5), see for example [20, Chap. II, Sec. 8.3, Prop. 9]. Henceforth, we will denote without ambiguity by $\phi_i$ the electric potential on the surface $\Sigma_i$.
This regularity of the electric potential and normal currents leads to the well-known fact that the electric potential at any point $x \in \mathbb{R}^3$ is given by :

$$ \sigma(x)\phi(x) = \sigma_0\phi(\xi)(x) - \sum_{i=1}^{n+1} \frac{\sigma_i^- - \sigma_i^+}{4\pi} \int_{\Sigma_i} \phi_i(y)\nu_i(y) \cdot \nabla_y \left( \frac{1}{|x-y|} \right) d\mathcal{H}_i(y) \quad (10) $$

where $\mathcal{H}_i$ is the 2-dimensional Hausdorff measure on the surface $\Sigma_i$, see for example [10]. Note that the above formula is valid for $x \in \Sigma_i$ by taking the non-tangential limit approaching from int $\Sigma_i$. Using (10) we have that on each interface $\Sigma_k, k = 1, 2, \dots, n+1$ the first regularity condition of (9) can be rewritten as

$$ \frac{\sigma_k^+ + \sigma_k^-}{2} \phi_k^{\pm}(x) = \sigma_0 \phi(\xi)(x) - \sum_{i=1, i \neq k}^{n+1} (\sigma_i^- - \sigma_i^+) K_i \phi_i(x) - (\sigma_k^- - \sigma_k^+) K_k \phi_k(x), \quad (11) $$

where $K_i$ is the double layer potential defined on the interface $\Sigma_i$. Once (11) is satisfied the continuity of the normal derivatives of the double layer potentials across the interfaces ensures that the second condition of (9) is also satisfied. Due to the condition that the electric conductivity outside $\Omega$ is zero, the electric potential has to satisfy,

$$ \sigma_n \partial_{\nu_{n+1}} \phi^{-}(x) = 0, \quad (12) $$

for $x \in \Sigma_{n+1}$. It is important to note that for one to be able to obtain the electric potential associated with $\xi$, it is sufficient to solve a system of equations with (11) on the surfaces $\Sigma_k, k = 1, 2, \dots, n$ and (12) on the surface $\Sigma_{n+1}$.

Note that the electric potentials on the surfaces $\Sigma_i$ also produce magnetic fields that have to be considered when looking at the magnetic flux density associated with the primary current $\xi$. At any point $x \in \mathbb{R}^3$ have the following expression for the magnetic flux density associated with the primary current $\xi$,

$$ \mathbf{B}(x) = \mathbf{B}(\xi)(x) - \mu \sum_{i=1}^{n+1} \frac{\sigma_i^- - \sigma_i^+}{4\pi} \int_{\Sigma_i} \nu_i(y) \times \nabla_y \frac{1}{|x-y|} \phi_i(y) d\mathcal{H}_i(y), \quad (13) $$

see for example [10, Eq. (17)], where the $\phi_i$'s on the surfaces, $\Sigma_i$, are the same as the surface potentials in (10). Note that in (13) care need to be taken when $x \in \Sigma_i$, see for example in the proof of Proposition 1 below.

## 4 [$W^{1/2,2}]^3$ sources

In what follows $\Sigma_0$ denotes a closed Lipschitz surface inside $\Omega_0$, which is the support of the sources. We now take $[\Xi]^3$ to be the subspace of $[W^{1/2,2}(\Sigma_0)]^3$ composed of vector-fields normally oriented to $\Sigma_0$, that is we take vector-fields $M_{\Sigma_0} \in [W^{1/2,2}(\Sigma_0)]^3$ of the form $M_{\Sigma_0} = M_{\Sigma_0} \nu_0$, with $M_{\Sigma_0} \in W^{1/2,2}(\Sigma_0)$. Both the orientation and magnitude of $M_{\Sigma_0}$ are encoded by $M_{\Sigma_0} \in W^{1/2,2}(\Sigma_0)$, we take $M_{\Sigma_0}$ to be oriented in the same
direction as $\nu_0$ when $M_{\Sigma_0}$ is positive and in the opposite direction otherwise. Note that (6) can be rewritten as

$$ 
\begin{aligned}
\sigma_0 \phi(\mathbf{M}_{\Sigma_0})(x) &= \frac{1}{4\pi} \int_{\Sigma_0} \mathbf{M}_{\Sigma_0}(y) \cdot \frac{(x-y)}{|x-y|^3} d\mathcal{H}_0(y), \\
&= \frac{1}{4\pi} \int_{\Sigma_0} M_{\Sigma_0}(y) \frac{(x-y)}{|x-y|^3} \cdot \nu_0(y) d\mathcal{H}_0(y), \quad (14) \\
&= K_0 M_{\Sigma_0}(x),
\end{aligned}
$$

for $x \in \mathbb{R}^3 \setminus \Sigma_0$. From [17, Thm 4.1] we have that $\phi \in W^{1,2}(\Omega)$ and $\phi \in W^{1,2}(\mathbb{R}^3 \setminus \Omega)$. For $x \in \Sigma_0$, we have that by approaching $x$ non-tangentially

$$ \sigma_0 \phi(\mathbf{M}_{\Sigma_0})(x) = \pm \frac{M_{\Sigma_0}(x)}{2} + K_0 M_{\Sigma_0}(x), $$

where the $-$ and $+$ are from approaching the boundary from interior and exterior of $\Sigma_0$, respectively.

## 4.1 Forward model for electric potential

We can now rewrite (10) as

$$ 
\begin{aligned}
\sigma(x)\phi(x) &= \sigma_0\phi(\mathbf{M}_{\Sigma_0})(x) - \sum_{i=1}^{n+1} \frac{\sigma_i^- - \sigma_i^+}{4\pi} \int_{\Sigma_i} \phi_i(y)\nu_i(y) \cdot \nabla_y \left( \frac{1}{|x-y|} \right) d\mathcal{H}_i(y) \quad (15) \\
&= K_0 M_{\Sigma_0}(x) - \sum_{i=1}^{n+1} (\sigma_i^- - \sigma_i^+) K_i \phi_i(x),
\end{aligned}
$$

We can therefore rewrite (11) as

$$ \frac{\sigma_k^+ + \sigma_k^-}{2} \phi_k^{\pm}(x) = K_0 M_{\Sigma_0}(x) - \sum_{i=1, i \neq k}^{n+1} (\sigma_i^- - \sigma_i^+) K_i \phi_i(x) - (\sigma_k^- - \sigma_k^+) K_k \phi_k(x). \quad (16) $$

We now discuss the numerical implementation of (15) and (16) on triangular meshes of the surfaces $\Sigma_i$. Since we have made the assumption that $\Omega$ is a Lipschitz domain then $\mathcal{K}: L^2(\Sigma) \to L^2(\Sigma)$ fails generally to be compact, see for example [15]. When $\mathcal{K}: L^2(\Sigma) \to L^2(\Sigma)$ is compact we have that $(\pm \frac{1}{2}I + \mathcal{K}): L^2(\Sigma) \to L^2(\Sigma)$ has a canonical representation, that is,

$$ (\pm \frac{1}{2}I + \mathcal{K})f = \sum_{j \ge 1} \alpha_j \langle f, u_j \rangle u_j, $$

where $\{\alpha_j\}$ and $\{u_j\}$ are the eigenvalues and an orthonormal basis of eigenfunctions of $(\pm \frac{1}{2}I + \mathcal{K})$, respectively. This is true for example when $\Sigma$ is a sphere in which case the $u_j$ are spherical harmonics, see for example [21], hence we can derive explicit expressions for (15). In general, if $\Sigma$ is $C^1$ then $\mathcal{K}$ is a compact operator, see for example [15],
hence a canonical representation as was given above is achievable. It then is possible to numerically approximate the eigenvalues and eigenfunctions in these smooth cases, see for example [22]. In these cases we can use the method of fundamental solutions to build (15) and (16), see [8] for an example of this application of the method.

For $\Sigma$ a Lipschitz surface we propose the following use of boundary elements methods for (15) and (16). Given a surface $\Sigma$, we begin by triangulating the surface to obtain $\Sigma_T$. On $\Sigma_T$ we will represent each function $f \in W^{1/2,2}(\Sigma)$ by considering its values on the vertices of $\Sigma_T$. We assume that on each triangle the function can be represented by linear shape functions, that is, on each triangle there are three linear basis functions each of which has value one on one vertex and zero on the other two. Thus, given a function $f \in W^{1/2,2}(\Sigma)$, on each triangle $T_k$ of $\Sigma_T$ we write,

$$f(y) = \sum_{j=1}^{3} f_{kj} \psi_{kj}(y),$$

where $y \in T_k$, $f_{kj}$ is the value of $f$ on the $j$-th vertex of the triangle $T_k$ and $\psi_{kj}$ is the linear shape function on $T_k$ that has value one on the $j$-th vertex of the triangle. Now given a point $x \in \mathbb{R}^3$ we wish to compute the double layer potential $Kf(x)$ for $f \in W^{1/2,2}(\Sigma)$. We use the analytic formulation proposed in [23], which enable us to numerically compute exactly the double layer potential defined on $\Sigma_T$ even when $x \in \Sigma_T$. With this formulation we write

$$Kf(x) = \mathbf{H}(x)\mathbf{f},$$

where $\mathbf{H}(x)$ is a row vector in which the $l$-th element of $\mathbf{H}(x)$ is the sum of the contribution to the double layer potential of each triangle that has the $l$-th vertex of the triangulation as a vertex and $\mathbf{f}$ is a column vector of the values of $f$ on the vertices of $\Sigma_T$ with the $l$-th element of $\mathbf{f}$ is the value of $f$ on the $l$-th vertex of the triangulation. Thus depending on where $x$ is located, $\mathbf{H}(x)\mathbf{f}$ is either $Kf(x)$ or $(-\frac{1}{2}I + K)f(x)$, here we have a “-” for $x \in \Sigma_T$ as we assume that the approach is from the interior. Using this notation, (15) can then be written as

$$\sigma(x)\phi(x) = \sigma_0\mathbf{H}_0(x)\mathbf{\Phi}_0 - \sum_{i=1}^{n+1}(\sigma_i^- - \sigma_i^+)\mathbf{H}_i(x)\mathbf{\Phi}_i,$$

where $\mathbf{H}_i(x)$ is as described above and also depends on the surface $\Sigma_{i,T}$, the triangulation of $\Sigma_i$ and $\mathbf{\Phi}_i$ are the values of $\phi_i$ at the vertices of $\Sigma_{i,T}$, respectively. For $x \in \Sigma_{k,T}$, $k = 1, 2, \dots, n+1$, (16) becomes

$$\sigma_k^-(x)\phi(x) + (\sigma_k^--\sigma_k^+)\mathbf{H}_k(x)\mathbf{\Phi}_k = \sigma_0\mathbf{H}_0(x)\mathbf{\Phi}_0 - \sum_{i=1, i \neq k}^{n+1}(\sigma_i^--\sigma_i^+)\mathbf{H}_i(x)\mathbf{\Phi}_i.$$

We can use the same idea as above to compute the gradient of the double layer potential at any point $x$ for any function $f \in W^{1/2,2}(\Sigma)$. To be able to solve the
forward cortical mapping problem we only need the normal derivatives of the double layer potential at the outer most surface. We will have the following formulation

$$ \partial_{\nu} K f(x) = \mathbf{N}(x) \mathbf{f}, $$

where $\mathbf{N}(x)$ is constructed in the same manner as $\mathbf{H}(x)$. It then follows that the normal derivative of (15) for $x \in \Sigma_{n+1,T}$ becomes

$$ \sigma_0 \mathbf{N}_0(x) \mathbf{\Phi}_0 = \sum_{i=1}^{n+1} (\sigma_i^- - \sigma_i^+) \mathbf{N}_i(x) \mathbf{\Phi}_i, $$

Now given the surfaces $\Sigma_{i,T}$, $\mathbf{\Phi}_0$ and letting $X_i$ the set of all vertices on all surfaces $\Sigma_{i,T}$ we can build the linear system

$$ \left\{ \left( \sigma_k^- \mathbf{I} + (\sigma_k^- - \sigma_k^+) \mathbf{H}_k(X_k) \right) \mathbf{\Phi}_k + \sum_{i=1, i \neq k}^{n+1} (\sigma_i^- - \sigma_i^+) \mathbf{H}_i(X_k) \mathbf{\Phi}_i = \sigma_0 \mathbf{H}_0(X_k) \mathbf{\Phi}_0 \right\}_{k=1,2,...,n} \\ \sum_{i=1}^{n+1} (\sigma_i^- - \sigma_i^+) \mathbf{N}_i(X_{n+1}) \mathbf{\Phi}_i = \sigma_0 \mathbf{N}_0(X_{n+1}) \mathbf{\Phi}_0. \tag{17} $$

Solving the above linear system for $\mathbf{\Phi}_k$, $k = 1, 2, ..., n + 1$, solves the forward cortical mapping problem.

As alluded to earlier to compute the matrices $\mathbf{H}$ and $\mathbf{N}$ we will implement a method suggested in [23] which results in exact expressions for the quantities that are required to build the matrices for the meshed surfaces, $\Sigma_{i,T}$. This allows to compute these matrices for any arbitrary point $x \in \mathbb{R}^3$ hence allowing the numerical approximations required for either EEG or sEEG with a high accuracy.

## 4.2 Forward model for magnetic flux density

We now turn our attention to the magnetic flux density associated with vector-fields in $[W^{1/2,2}(\Sigma_0)]^3$ that were introduced in the previous subsection. It follows from (7) that the magnetic vector potential of $\mathbf{M}_{\Sigma_0}$ is given by

$$ \mathbf{A}(\mathbf{M}_{\Sigma_0})(x) = \frac{\mu}{4\pi} \int_{\Sigma_0} \frac{\mathbf{M}_{\Sigma_0}(y)}{|x-y|} d\mathcal{H}(y) = \mu S \mathbf{M}_{\Sigma_0}(x), \quad x \in \mathbb{R}^3 \setminus \Sigma_0, \tag{18} $$

with the integral being taken in the principal value sense for $x \in \Sigma_0$, see for example [19, Eq. (5.32)]. Note that the last equality of (18) is a direct application of the single layer potential hence $\mathbf{A}(\mathbf{M}_{\Sigma_0})$ has components that are equal to the single layer potentials of the corresponding components of $\mathbf{M}_{\Sigma_0}$. Since $W^{1/2,2}(\Sigma_0) \subset L^2(\Sigma_0) \subset W^{-1/2,2}(\Sigma_0)$ we have from [17, Thm 3.1] $\mathbf{A}(\mathbf{M}_{\Sigma_0}) \in [W^{1,2}(\Omega)]^3$ and $\mathbf{A}(\mathbf{M}_{\Sigma_0}) \in [W^{1,2}(\mathbb{R}^3 \setminus \bar{\Omega})]^3$.

Since $\mathbf{B} = -\nabla \times \mathbf{A}$ it follows that $\mathbf{B}(\mathbf{M}_{\Sigma_0}) \in [L^2(\Omega)]^3$ and $\mathbf{B}(\mathbf{M}_{\Sigma_0}) \in [L^2(\mathbb{R}^3 \setminus \bar{\Omega})]^3$ and the magnetic flux density associated with a vector-field $\mathbf{M}_{\Sigma_0} \in [W^{1/2,2}(\Sigma_0)]^3$ is given by

$$ \mathbf{B}(\mathbf{M}_{\Sigma_0})(x) = \frac{\mu}{4\pi} \int_{\Sigma_0} \mathbf{M}_{\Sigma_0}(y) \times \nabla_y \frac{1}{|x-y|} d\mathcal{H}(y) = \mu \nabla \times S \mathbf{M}_{\Sigma_0}(x), \quad x \in \mathbb{R}^3 \setminus \Sigma_0. $$
Since $\mathbf{M}_{\Sigma_0} = M_{\Sigma_0} \nu_0$ we can rewrite the above equation as

$$ \mathbf{B}(\mathbf{M}_{\Sigma_0})(x) = \frac{\mu}{4\pi} \int_{\Sigma_0} \nu_0(y) \times \nabla_y \frac{1}{|x-y|} M_{\Sigma_0}(y) d\mathcal{H}(y) = \mu \mathbf{S} M_{\Sigma_0}(x), \quad x \in \mathbb{R}^3 \setminus \Sigma_0. \qquad (19) $$

Note that the last equality of (19) is how we define the operator $\mathbf{S}: W^{1/2,2}(\Sigma_0) \to \mathbb{R}^3$, which will be carried over to the discrete case. It follows from (13) that

$$ \begin{aligned} \mathbf{B}(x) &= \mathbf{B}(\mathbf{M}_{\Sigma_0})(x) - \mu \sum_{i=1}^{n+1} \frac{\sigma_i^- - \sigma_i^+}{4\pi} \int_{\Sigma_i} \nu_i(y) \times \nabla_y \frac{1}{|x-y|} \phi_i(y) d\mathcal{H}_i(y) \\ &= \mathbf{S} M_{\Sigma_0}(x) - \mu \sum_{i=1}^{n+1} \frac{\sigma_i^- - \sigma_i^+}{4\pi} \int_{\Sigma_i} \nu_i(y) \times \nabla_y \frac{1}{|x-y|} \phi_i(y) d\mathcal{H}_i(y), \end{aligned} \qquad (20) $$

where the $\phi_i$'s on the surfaces are surface potentials obtained from the forward model of the electric potential associated with $M_{\Sigma_0}$. Note that MEG measures $\mathbf{B} \cdot \mathbf{v}$, for some known vector-field $\mathbf{v}$. Typically, $\mathbf{v}$ is taken to be a radial vector-field, hence if all of the $\Sigma_i$ above are spherical, all the vector-field of the form $M_{\Sigma_0} \nu_0$ will result in null MEG measurements, see [10, Eq. (20)], which make a spherical head model uninteresting for MEG in this context.

We now look at how to numerically compute (20) for Lipschitz surfaces $\Sigma_i$. To that end we need to be able to express $\mathbf{S} M_{\Sigma_0}$ as a linear combination of the $\mathbf{S} M_{\Sigma_{0i}}$'s where the $M_{\Sigma_{0i}} \in W^{1/2,2}(\Sigma_0)$ are elements of the basis of $W^{1/2,2}(\Sigma_0)$. Note that it only suffices that we look at how to numerically compute the expression given in (19) because to compute (20) we repeatedly apply the same idea. We begin with a triangulation of the surfaces and discretisation of $M_{\Sigma} \in W^{1/2,2}(\Sigma)$ as was done in the previous subsection. Note that (19) on $\Sigma_T$ is given by

$$ \mathbf{B}(x) = \sum_{T_k} \frac{\mu}{4\pi} \int_{T_k} \nu_{T_k}(y) \times \nabla_y \frac{1}{|x-y|} M_{\Sigma}(y) d\mathcal{H}(y), $$

since $\nu_{T_k}$ is taken to be constant on each triangle hence we have that

$$ \begin{aligned} \mathbf{B}(x) &= \mu \sum_{T_k} \nu_{T_k} \times \frac{1}{4\pi} \int_{T_k} \nabla_y \frac{1}{|x-y|} M_{\Sigma}(y) d\mathcal{H}(y) \\ &= \mathbf{S}(x)\mathbf{M}_{\Sigma}, \end{aligned} $$

where $\mathbf{S}(x)$ is a row vector that depends on $x$ which is generated in a similar manner as $\mathbf{H}(x)$ and $\mathbf{M}_{\Sigma}$ is a column vector of the values of $M_{\Sigma}$ on the vertices of $\Sigma_T$, where we use the formula given in [23] to compute the integral above on each triangle. Hence (20) is written in discrete form as

$$ \mathbf{B}(x) = \mathbf{S}_0(x)\mathbf{\Phi}_0 - \sum_{i=1}^{n+1} (\sigma_i^- - \sigma_i^+) \mathbf{S}_i(x)\mathbf{\Phi}_i, $$

where the $\mathbf{\Phi}_i$, $i = 1, 2, \dots, n+1$, are obtained from the forward model for electrical potential.
## 4.3 Silent Sources

We will conclude this section by discussing silent sources. We begin with a definition:

**Definition 1.** Let $E \subset \mathbb{R}^3\text{supp}(\xi)$ then a non-zero vector-field $\xi \in [\Xi]^3$ is electrically (magnetically) silent in $E$ if it produces an identically zero electric potential (magnetic flux density) in $E$.

The existence of silent sources guarantees the non-uniqueness of solutions to the source identification problems we aim to solve. It is therefore important to understand the nature of the silent sources as this may help mitigate their impact on the uniqueness of solutions.

As an example, take $[\Xi]^3$ to be the subspace of $[W^{1/2,2}(\Sigma_0)]^3$ composed of vector-fields of the form $M_{\Sigma_0}\nu_0$, with $M_{\Sigma_0} \in W^{1/2,2}(\Sigma_0)$. We will now explore the silent sources among these vector-fields.

**Proposition 1.** Vector-fields of the form $M_{\Sigma_0}\nu_0$, with $M_{\Sigma_0} \in W^{1/2,2}(\Sigma_0)$ such that $M_{\Sigma_0}$ is constant are both electrically and magnetically silent in $\mathbb{R}^3\setminus\Omega_0$.

*Proof.* From (14) we have that

$$ \sigma_0\phi(\mathbf{M}_{\Sigma_0})(x) = K_0M_{\Sigma_0}(x), \quad x \in \mathbb{R}^3\setminus\Omega_0, $$

and from [16, Ex. 6.14] we have that the above equation is identically zero in $\mathbb{R}^3\setminus\Omega_0$ if $M_{\Sigma_0}$ is a constant. Hence the proposition is proved for the electric silence.

It also follows from [24, Lem. 4.3] that

$$ \mathbf{B}(\mathbf{M}_{\Sigma_0})(x) = \mathbf{S}M_{\Sigma_0}(x) = \nabla \times S(M_{\Sigma_0}\nu)(x) = S(\nu \times \nabla M_{\Sigma_0})(x), $$

where $\nabla M_{\Sigma_0}$ is the gradient of $M_{\Sigma_0}$ on $\Sigma_0$. It follows from [15, Thm 3.3] that all $\mathcal{S}: L^2(\Sigma_0) \to L^2(\Sigma_0)$ is injective. Hence, if $\mathcal{S}(M)$ is identically zero in $\mathbb{R}^3\setminus\Omega_0$ then its non-tangential limit to $\Sigma_0$ is identically zero on $\Sigma_0$ thus $M = 0$ on $\Sigma_0$. We can therefore conclude that, if $M_{\Sigma_0} \in L^2(\Sigma_0)$ is such that $\nabla M_{\Sigma_0}$ vanishes on $\Sigma_0$ then the resulting magnetic flux density is identically zero. Hence, all $M_{\Sigma_0}\nu_0$ such that $M_{\Sigma_0}$ is constant are magnetically silent in $\mathbb{R}^3\setminus\Omega_0$. $\square$

**Remark 1.** From [16, Ex. 6.14] we observe that if $M_{\Sigma_0} \in W^{1/2,2}(\Sigma_0)$ is constant then

$$ \sigma_0\phi(\mathbf{M}_{\Sigma_0})(x) = K_0M_{\Sigma_0}(x) = M_{\Sigma_0}, \quad x \in \text{int } \Sigma_0, $$

hence these vector-fields, $M_{\Sigma_0}\nu_0$, are not silent in $\Omega_0$. We can exploit this fact by combining EEG or MEG data with sEEG data. The sEEG data will in principle allow to eliminate the silent sources outlined in Proposition 1 in the source recovery.

As a second example take a closed set $S \subset \text{int } \Sigma_0$ and let $S$ be a slender set, that is, $\mathcal{L}_3(S) = 0$, where $\mathcal{L}_3$ is the Lebesgue measure on $\mathbb{R}^3$ and $\mathcal{L}_3(\mathbb{R}^3\setminus S) = \infty$. We will take $[\Xi]^3 = [\mathcal{M}(S)]^3 \subset [\mathcal{M}(\text{int } \Sigma_0)]^3$, that is, vector-valued measures supported on the
slender set $S$. For $\xi \in [\mathcal{M}(\text{int } \Sigma_0)]^3$ we define $\nabla \cdot \xi$ in the sense of distributions and when $\nabla \cdot \xi = 0$ (divergence-free) we mean

$$ \int_{\text{int } \Sigma_0} \nabla u \cdot d\xi = \int_{\text{int } \Sigma_0} \sum_{j=1}^{3} \frac{\partial u}{\partial x_j} d\xi_j = 0, $$

for all $u \in C_c^\infty(\text{int } \Sigma_0)$.

**Proposition 2.** Every $\xi \in [\mathcal{M}(\text{int } \Sigma_0)]^3$ that has a slender support and is divergence-free is electrically silent in $\mathbb{R}^3 \setminus \Omega_0$.

*Proof.* This is a direct application of [25, Thm 2.2]. □

Examples of divergence-free vector-fields on slender sets in $\mathbb{R}^3$ can be constructed as follows. Let $\gamma : [0, l] \to \mathbb{R}^3$ be a Lipschitz mapping and let $S := \gamma([0, l])$. If $\gamma$ is such that

$$ \mathcal{H}_1(\gamma([a, b])) = b - a, \quad \forall [a, b] \subset [0, l], $$

then $\gamma$ is an orientable rectifiable curve. Note that $S$ is slender and on $S$ define the vector measure $\mathbf{R}_\gamma$ through the relation

$$ \langle \mathbf{R}_\gamma, \mathbf{f} \rangle = \int_0^l \mathbf{f}(\gamma(t)) \cdot \gamma'(t) dt, \quad \text{for } \mathbf{f} \in [C_c(\mathbb{R}^3)]^3, $$

where $\gamma'$ is the unit tangent vector of $\gamma$. Let the endpoints of $S$ be $s, t \in \mathbb{R}^3$, it follows that $\mathbf{R}_\gamma = \gamma' \mathcal{H}_1$ and the divergence of $\mathbf{R}_\gamma$ is given as

$$ \int_0^l \nabla u(\gamma(t)) \cdot \gamma'(t) dt = \int_\gamma \dot{u} d\mathcal{H}_1 = u(s) - u(t), \quad \forall u \in C_c^\infty(\text{int } \Sigma_0), $$

where $\dot{u}$ is the derivative of $u$ along $\gamma$. Hence if $s = t$ then $\mathbf{R}_\gamma$ is divergence free and by the proposition above is electrically silent. Concretely, take $S$ to be a circle embedded on $\mathbb{R}^3$ by taking

$$ \gamma := \begin{pmatrix} \sin t \\ \cos t \\ 0 \end{pmatrix}, $$

with $t \in [0, 2\pi]$. More details and example can be found in [26].

# 5 Inverse problems

In this section we consider the general situation when elements of $[\Xi]^3$ are vector-fields with components that are elements of a Banach space, $\Xi$, supported in $\bar{\Omega}_0$. When solving the inverse problems for MEG, EEG and sEEG, we need that the electric potential, $\phi$, associated with the recovered source satisfies the conditions (9) and (12). We recall (11) and (12) here

$$ -\sigma_0 \phi(\xi)(x) + \frac{\sigma_k^+ + \sigma_k^-}{2} \phi_k + \sum_{i=1, i \neq k}^{n+1} (\sigma_i^- - \sigma_i^+) K_i \phi_i(x) - (\sigma_k^- - \sigma_k^+) K_k \phi_k(x) = 0, \quad (21) $$

$$ \sigma_n \partial_{\nu_{n+1}} \phi(x) = 0, $$
and let $\mathcal{C} : [\Xi]^3 \times L^2(\Sigma_1) \times \dots \times L^2(\Sigma_{n+1}) \to L^2(\Sigma_1) \times L^2(\Sigma_2) \times \dots \times L^2(\Sigma_n) \times L^2(\Sigma_{n+1})$ be the LHS of (21). Note that the null space of $\mathcal{C}$, which we denote $C_0$, consists of those surface electric potentials and their associated sources that satisfy the conditions (9) and (12). In other words, we look for solutions to inverse problem of MEG, EEG and/or sEEG in $C_0$. For the discrete setting we rewrite (17) as

$$ \left\{ \begin{aligned} & -\sigma_0 \mathbf{\Phi}(\boldsymbol{\xi})(X_k) + (\sigma_k^- \mathbf{I} + (\sigma_k^- - \sigma_k^+) \mathbf{H}_k(X_k)) \mathbf{\Phi}_k \\ & \qquad + \sum_{i=1, i \neq k}^{n+1} (\sigma_i^- - \sigma_i^+) \mathbf{H}_i(X_k) \mathbf{\Phi}_i = \mathbf{0} \\ & -\sigma_0 \partial_{\nu_{n+1}} \mathbf{\Phi}(\boldsymbol{\xi})(X_{n+1}) + \sum_{i=1}^{n+1} (\sigma_i^- - \sigma_i^+) \mathbf{N}_i(X_{n+1}) \mathbf{\Phi}_i = \mathbf{0}, \end{aligned} \right. \quad (22) $$

where $\sigma_0 \mathbf{\Phi}(\boldsymbol{\xi})(X_k)$ and $\sigma_0 \partial_{\nu_{n+1}} \mathbf{\Phi}(\boldsymbol{\xi})(X_{n+1})$ are appropriate discretisation of (6) and $\sigma_0 \partial_{\nu_{n+1}} \phi(x)$, respectively. Hence, (22) gives the discrete version of $\mathcal{C}$ on the discretised surfaces $\Sigma_{i,T}$, $i = 1, 2, \dots, n+1$. The LHS of (22) will be denoted $\mathbf{C}$ and the null space of $\mathbf{C}$ by $C_0$.

## 5.1 sEEG, EEG and MEG problems

In the inverse source localisation problems, the electric potential ($\phi$) or magnetic flux density ($\mathbf{B}$) is used. $\phi$ is measured intra-cranially (sEEG) or on the scalp (EEG) whereas $\mathbf{B}$ is measured close to the scalp (OPM-MEG) or at a distance away from the head (SQUID-MEG). In any instance, we wish to solve the problem that given point-wise measurements of the electric potential ($\phi$) or magnetic flux density ($\mathbf{B}$) in some subset of $\mathbb{R}^3$, find $\boldsymbol{\xi} \in [\Xi]^3$ and $\phi_i \in L^2(\Sigma_i)$, $i = 1, 2, \dots, n+1$, such that (10) or (20), respectively, are satisfied.

Due to the existence of silent sources for sEEG, EEG and MEG as highlighted in Section 4.3, we can only recover $\boldsymbol{\xi} \in [\Xi]^3$ up to silent sources. Further, since we only have point-wise data for sEEG, EEG and MEG, the unique recovery of $\boldsymbol{\xi} \in [\Xi]^3$ is also negatively impacted by the existence of $\boldsymbol{\xi} \in [\Xi]^3$ that generate an electric potential or magnetic flux density that vanishes at the measurement points. In practice, $\mathbf{B}$ is known only point-wise and only a component of it is considered with the additional issue that the electric potential $\phi$ in $\bar{\Omega}$ is also unknown. Due to the non-uniqueness of solutions highlighted here these inverse problems are ill-posed, hence we solve Tikhonov regularised problems.

Keeping in mind that the electric potential that is reconstructed from the point-wise measurements has to have the regularity stipulated in (9) and (12). As indicated earlier we look for solutions in the null space $C_0$, which is exactly solving the cortical mapping problem. We build the problems in such a way that they incorporate the regularity requirements of the potential and normal current. There are a multitude of ways this can be achieved, for example, in [4] the authors used a projector onto $C_0$ and looked for solutions directly in $C_0$. In our case we will use the projection onto $C_0$ as a regulariser in the associated Tikhonov problems.
## 5.2 Existence of solutions to the sEEG, EEG and MEG problems

As highlighted above when solving the inverse problems for sEEG, EEG and MEG we aim to recover the source $\xi \in [\Xi]^3$ and by extension the surface electric potentials $\phi \in L^2(\Sigma_i)$, $i = 1, 2, \dots, n+1$. We now mathematically set up these problems and show that the solutions exist. To that end let $\mathcal{S}$ be the product Banach space

$$[\Xi]^3 \times L^2(\Sigma_1) \times L^2(\Sigma_2) \times \dots \times L^2(\Sigma_{n+1})$$

endowed with the norm

$$\| \cdot \|_{[\Xi]^3} + \| \cdot \|_{L^2(\Sigma_1)} + \| \cdot \|_{L^2(\Sigma_2)} + \dots + \| \cdot \|_{L^2(\Sigma_{n+1})},$$

and let

$$\mathcal{D}_1 = L^2(D),$$

be the Hilbert space which corresponds to the data measured with EEG, MEG and sEEG measurements with $D \subset \mathbb{R}^3$ depending on the measurement modalities used. Finally we would let the product Hilbert space

$$\mathcal{D}_2 = L^2(\Sigma_1) \times L^2(\Sigma_2) \times \dots \times L^2(\Sigma_n) \times L^2(\Sigma_{n+1}),$$

endowed with the norm

$$(\| \cdot \|_{L^2(\Sigma_1)}^2 + \| \cdot \|_{L^2(\Sigma_2)}^2 + \dots + \| \cdot \|_{L^2(\Sigma_{n+1})}^2)^{1/2},$$

which is useful in the study of (21). Let $\mathcal{F}_i: \mathcal{S} \to \mathcal{D}_i$ for $i = 1, 2$ be linear operators with $\mathcal{F}_1$ as the forward model of EEG, MEG and/or sEEG, that is, the formulas given in by (10), (13) and/or (10), respectively and $\mathcal{F}_2 = \mathcal{C}$. Thus, when solving the inverse problems for MEG, EEG and/or sEEG we consider the appropriate functional

$$\mathcal{T}_{f,\alpha,\beta,\lambda_0,\dots,\lambda_{n+1}}(\xi, \phi_1, \phi_2, \dots, \phi_{n+1}) := \alpha \|\mathcal{F}_1(\xi, \phi_1, \phi_2, \dots, \phi_{n+1}) - f\|_{\mathcal{D}_1}^2 \\ + \beta \|\mathcal{F}_2(\xi, \phi_1, \phi_2, \dots, \phi_{n+1})\|_{\mathcal{D}_2}^2 + \lambda_0 R(\|\xi\|_{[\Xi]^3}) + \sum_{j=1}^{n+1} \lambda_j \|\phi_j\|_{L^2(\Sigma_j)}^2, \quad (23)$$

with $f$ measured data and $R: [0, \infty) \to [0, \infty)$ a convex function. Thus we solve the following problem:

**Problem 1.** Given data $f \in \mathcal{D}_1$ and $\alpha, \beta, \lambda_j > 0$ find $(\xi, \phi_1, \phi_2, \dots, \phi_{n+1})_\lambda \in \mathcal{S}$ such that

$$(\xi, \phi_1, \phi_2, \dots, \phi_{n+1})_\lambda = \underset{(\xi, \phi_1, \phi_2, \dots, \phi_{n+1}) \in \mathcal{S}}{\arg\inf} \mathcal{T}_{f,\alpha,\beta,\lambda_0,\dots,\lambda_{n+1}}(\xi, \phi_1, \phi_2, \dots, \phi_{n+1}).$$

**Remark 2.** We will prove the existence of a solution to Problem 1 in Theorem 1 below keeping in mind that we are mostly interested in $[\Xi]^3$ that is either $[W^{1/2,2}(\Sigma_0)]^3$ or $[\mathcal{M}(\Sigma_0)]^3$ or $[\mathcal{M}(\text{int} \Sigma_0)]^3$, with $\mathcal{M}$ being the Banach space of measures endowed with the total variation norm. By using embeddings of various Sobolev space into the
space of continuous function we view spaces of measures as being contained in the duals of certain Sobolev spaces. We take the view that $[M(\text{int } \Sigma_0)]^3 \subset [(W^{1,q}(\text{int } \Sigma_0))^*]^3$ for $q > 3$ from the Sobolev embedding theorem, see for example [27, Thm 5.4, Part II] hence we will discuss about the Newton potential in these spaces. We also have $[M(\Sigma_0)]^3 \subset [(W^{1-1/q,q}(\Sigma_0))^*]^3$ for appropriate choices of $q$, see for example [28, Thm 4.57, Thm 4.58] states that $W^{1,q}(\Omega)$ embeds into a space of functions continuous on $\bar{\Omega}$ hence the traces $W^{1-1/q,q}(\Sigma_0)$ are continuous. In the case of $[W^{1/2,q}(\Sigma_0)]^3$ and $[M(\Sigma_0)]^3$ we only need to look at how the layer potentials behave.

**Theorem 1.** A unique solution to Problem 1 exists.

*Proof.* We use the result [29, Thm 3.1] to make this conclusion hence we need only show that Problem 1 satisfies the assumptions of [29, Thm 3.1].

(IP1) We note that the duals of $\mathcal{S}$ and $\mathcal{D}_2$ are isometric to $([\Xi]^3)^* \times L^2(\Sigma_1) \times L^2(\Sigma_2) \times \cdots \times L^2(\Sigma_{n+1})$ and $L^2(\Sigma_1) \times L^2(\Sigma_2) \times \cdots \times L^2(\Sigma_{n+1})$, respectively, and we will associate $\mathcal{S}$ with its weak* topology, $\mathcal{D}_1$ and $\mathcal{D}_2$ with their weak topologies.

(IP2) From [30, Prop. 3.5] we have that norms are weakly lower semi-continuous hence the norm of $\mathcal{D}_1$ and $\mathcal{D}_2$ are weakly lower semi-continuous.

(IP3) We will discuss the continuity of $\mathcal{F}_1$ and $\mathcal{F}_2$.

(i) We discuss the continuity of $\mathcal{F}_1$ corresponding to (10) and $\mathcal{F}_2$ which corresponds to (21) because of the similarities.

(a) When $\Xi = W^{s,p}(\partial\Omega)$ for $0 < s < 1$, $1 < p < \infty$, the double layer potential from $W^{s,p}(\partial\Omega)$ into $W^{1+1/p-s,p}(\Omega)$ for $0 < s < 1$, $1 < p < \infty$, is continuous, see [17, Thm 4.1].

(b) For the case $\Xi = (W^{1-s,p}(\partial\Omega))^*$, we begin by noting that (6) can be rewritten

$$ \sigma\phi(\xi)(x) = \frac{1}{4\pi} \int_{\partial\Omega} \frac{\xi(y) \cdot (x - y)}{|x - y|^3} d\mathcal{H}(y), $$

which can be seen as the sum of the quantities

$$ \frac{1}{4\pi} \int_{\partial\Omega} \frac{\xi_i(y)(x - y)_i}{|x - y|^3} d\mathcal{H}(y), $$

for $i = 1, 2, 3$, each of which is the $i$-th component of the gradient of the single layer potential $S\xi_i(x)$. The single layer potential maps $W^{-s,p}(\partial\Omega)$ into $W^{1+1/p-s,p}(\Omega)$ for $0 < s < 1$, $1 < p < \infty$ continuously and $\nabla$ maps $W^{s,p}(\Omega)$ to $W^{s-1,p}(\Omega)$ continuously for $s > 0$, $1 < p < \infty$. Hence $\phi$ maps $(W^{1-s,p}(\partial\Omega))^*$ to $W^{1/q-s,q}(\Omega)$ for $0 < s < 1$, $1 < p, q < \infty$, $1/p+1/q = 1$, continuously.

(c) When $\Xi = (W^{1+s,q}(\Omega))^*$ for $-1 \le s \le 2$, $1 < p, q < \infty$, $1/p + 1/q = 1$, from the continuity of the Newton potential from the distributions $(W^{1+s,q}(\Omega))^*$ to $W^{1-s,p}(\Omega)$, see [17, Prop. 2.1], and the continuity of the double layer potential from $W^{s,p}(\partial\Omega)$ into $W^{1+1/p-s,p}(\Omega)$ with $\Omega$ a bounded Lipschitz domain.
(ii) We now discuss the continuity of $\mathcal{F}_1$ corresponding to (13).

(a) When $\Xi = W^{s,p}(\partial\Omega)$ for $0 \le s \le 1$ and $1 < p < \infty$, we begin by noting that

$$ \{(\nu \times \nabla)f : f \in W^{1-s,p}(\partial\Omega)\} \subseteq W^{-s,p}(\partial\Omega), $$

see [24, (4.7)]. From [24, Lem. 4.3] and its proof we conclude that each component of the terms that appear in the sum of (13) is continuous since [17, Thm 3.1] shows that the single layer potential is continuous from $W^{-s,p}(\partial\Omega)$ into $W^{1+1/p-s,p}(\Omega)$ for $0 < s < 1$, $1 < p < \infty$.

(b) When $\Xi = (W^{1-s,p}(\partial\Omega))^*$, we note that (7) can be rewritten as

$$ \mathbf{A}(\xi)(x) = \frac{\mu}{4\pi} \int_{\partial\Omega} \frac{\xi(y)}{|x-y|} d\mathcal{H}(y), $$

thus each component is a single layer potential. Since the single layer potential maps $W^{-s,p}(\partial\Omega)$ into $W^{1+1/p-s,p}(\Omega)$ for $0 < s < 1$, $1 < p < \infty$ continuously, we have that $\mathbf{B}$ maps $(W^{1-s,p}(\partial\Omega))^*$ to $W^{1/q-s,q}(\Omega)$ for $0 < s < 1$, $1 < p, q < \infty$, $1/p+1/q = 1$, continuously in each component.

(c) For the case $\Xi = (W^{1+s,q}(\Omega))^*$ for $-1 \le s \le 2$, $1 < p, q < \infty$, $1/p + 1/q = 1$, note that the first term on the RHS of (13) is the curl of the Newton potential. From (7) and (8) it follows that $\mathbf{B}$ is continuous from $(W^{1+s,q}(\Omega))^*$ to $W^{-s,p}(\Omega)$ in each component for $-1 \le s \le 2$, $1 < p, q < \infty$, $1/p + 1/q = 1$.

By choosing $s$ and $p$ appropriately we can conclude that $\mathcal{F}_1$ and $\mathcal{F}_2$, are continuous from $\mathcal{S}$ to $\mathcal{D}_1$ and $\mathcal{D}_2$, respectively in their weak topologies precised in (IP1).

(IP4) Note that

$$ \mathfrak{R} = \lambda_0 R(\|\xi\|_{[\Xi]^3}) + \sum_{j=1}^{n+1} \lambda_j \|\phi_j\|_{L^2(\Sigma_j)}^2 $$

is weak* lower semi-continuous in $\mathcal{S}$ from [30, Prop. 3.13]. Note that the zero element, $0$, of $\mathcal{S}$ is in the combined domain of $\mathcal{F}_1$ and $\mathcal{F}_2$, dom($\mathcal{F}_1, \mathcal{F}_2$), we have that $\mathfrak{R}(0) = 0$ hence $\mathfrak{R}$ is finite for at least one element of dom($\mathcal{F}_1, \mathcal{F}_2$), that is, $\mathfrak{R}$ is proper. Finally, $\mathfrak{R}$ is a convex combination of norms.

(IP5) We now show the that dom($\mathcal{F}_1, \mathcal{F}_2$) is weak* closed. We will again do this in two steps, first for the electrical potential and then for the magnetic flux density.

(i) Note that for the electric potential, the domain $[\Xi]^3$ is $[(W^{1+s,q}(\Omega))^*]^3$ or $[(W^{-s-1/p,q}(\partial\Omega))^*]^3$ for $-1 \le s \le 2$, $1 < p, q < \infty$ and these are weak* closed, see for example [30, Thm 3.33].

(ii) For the magnetic flux density we require that $[\Xi]^3$ is the set of elements $\xi \in [(W^{1+s,q}(\Omega))^*]^3$ for $-1 \le s \le 2$, $1 < p, q < \infty$ with $\nabla \times \xi \in [(W^{1+s,q}(\partial\Omega))^*]^3$ which is a closed subspace of $[(W^{1+s,q}(\Omega))^*]^3$ or $[\Xi]^3$ is $\xi \in [(W^{-s-1/p,q}(\partial\Omega))^*]^3$. In either case we have weak* closed sets, see for example [30, Thm 3.33].
Hence dom($\mathcal{F}_1, \mathcal{F}_2$) $\subset \mathcal{S}$ is weak* closed.

(IP6) Let $\gamma > 0$ and $\kappa > 0$ and set one of $\alpha, \beta, \lambda_0, \dots, \lambda_{n+1}$ equal to $\gamma$ define the set

$$ \mathcal{D} = \{ \mathbf{u} \in \mathrm{dom}(\mathcal{F}_1, \mathcal{F}_2) \subset \mathcal{S} : \mathcal{T}_{f,\alpha,\beta,\lambda_0,\dots,\lambda_{n+1}}(\mathbf{u}) \le \kappa \}. $$

Since $\mathcal{T}_{f,\alpha,\beta,\lambda_0,\dots,\lambda_{n+1}}(\mathbf{u})$ is continuous on $\mathrm{dom}(\mathcal{F}_1, \mathcal{F}_2) \subset \mathcal{S}$ it follows that $\mathcal{D}$ is a bounded set hence weak* compact from Banach-Alaoglu-Bourbaki Theorem, see for example [30, Thm 3.16].

Conditions (IP1)-(IP6) above show that Problem 1 satisfies the assumption of [29, Thm 3.1] hence a solution exists and the solution is unique due to the convexity of $\mathcal{T}_{f,\alpha,\beta,\lambda_0,\dots,\lambda_{n+1}}$. $\square$

**Remark 3.** Note that Problem 1 does not discuss the contamination of the data with noise. In practice $f \in \mathcal{D}_1$ is contaminated with noise and we assume that the noisy data $f^\delta \in \mathcal{D}_1$ and $\|f - f^\delta\|_{\mathcal{D}_1} \le \delta, \delta > 0$. The stability of Problem 1 and convergence of Problem 1 with respect to noise and regularisers $\alpha, \beta, \lambda_0, \dots, \lambda_{n+1}$ is also discussed in [29].

## 5.3 Alternating minimisation algorithm for solving inverse problems

We wish to solve Problem 1 by splitting it into a problem that solves for the source $\xi \in [\Xi]^3$ and a problem that solves for the surface electric potentials $\phi \in L^2(\Sigma_i)$. We iteratively solves these problems to obtain the solution we desire as discussed in the introduction of [31]. The alternating minimisation procedure is as follows, beginning with some initial guess

$$ (\xi^{\{0\}}, \phi_1^{\{0\}}, \phi_2^{\{0\}}, \dots, \phi_{n+1}^{\{0\}}), $$

then generate a sequence of solutions

$$ \{(\xi^{\{l\}}, \phi_1^{\{l\}}, \phi_2^{\{l\}}, \dots, \phi_{n+1}^{\{l\}})_\lambda\}_{l \in \mathbb{N}}, $$

by solving the following problems

$$
\begin{align}
\xi_{\lambda}^{\{l+1\}} &= \underset{\xi \in [\Xi]^3}{\arg\inf} \mathcal{T}_{f,\alpha,\beta,\lambda_0,\dots,\lambda_{n+1}}(\xi, \phi_1^{\{l\}}, \phi_2^{\{l\}}, \dots, \phi_{n+1}^{\{l\}}) \nonumber \\
(\phi_1^{\{l+1\}}, \phi_2^{\{l+1\}}, \dots, \phi_{n+1}^{\{l+1\}})_\lambda &\qquad (24) \\
&= \underset{(\phi_1, \dots, \phi_{n+1}) \in \mathcal{D}_2}{\arg\inf} \mathcal{T}_{f,\alpha,\beta,\lambda_0,\dots,\lambda_{n+1}}(\xi^{\{l\}}, \phi_1, \phi_2, \dots, \phi_{n+1}). \nonumber
\end{align}
$$

A closer inspection of $\mathcal{F}_1$ and $\mathcal{F}_2$ as we proposed be taken reveals that

$$
\begin{align*}
\mathcal{F}_1(\xi, \phi_1, \phi_2, \dots, \phi_{n+1}) &= \mathcal{F}_{1,1}(\xi) + \mathcal{F}_{1,2}(\phi_1, \phi_2, \dots, \phi_{n+1}) \\
\mathcal{F}_2(\xi, \phi_1, \phi_2, \dots, \phi_{n+1}) &= \mathcal{F}_{2,1}(\xi) + \mathcal{F}_{2,2}(\phi_1, \phi_2, \dots, \phi_{n+1}),
\end{align*}
$$
hence the problems that we solve in (24) have the variables $\xi \in [\Xi]^3$ and $(\phi_1, \phi_2, \dots, \phi_{n+1}) \in \mathcal{D}_2$ well separated in $\mathcal{F}_1$ and $\mathcal{F}_2$ hence we can implement methods that best recover each variable.

We now show that the sequence

$$ \{(\xi^{\{l\}}, \phi_1^{\{l\}}, \phi_2^{\{l\}}, \dots, \phi_{n+1}^{\{l\}})_{\lambda} \}_{l \in \mathbb{N}}, $$

generated by the above alternating minimisation algorithm results in the objective function converging linearly to the minimum of Problem 1. The result below is valid for all $[\Xi]^3$ outline in Remark 2 and Theorem 1.

**Theorem 2.** The sequence generated by the alternating minimisation algorithm converges to the minimum of the Tikhonov functional, $T_{f,\alpha,\beta,\lambda_0,\dots,\lambda_{n+1}}$, linearly.

*Proof.* Take

$$ \mathfrak{B}_1 = [\Xi]^3, $$

with

$$ \| \cdot \|_1 = \| \cdot \|_{[\Xi]^3}, $$

and

$$ \mathfrak{B}_2 = \mathcal{D}_2, $$

with

$$ \| \cdot \|_2 = (\| \cdot \|_{L^2(\Sigma_1)}^2 + \| \cdot \|_{L^2(\Sigma_2)}^2 + \dots + \| \cdot \|_{L^2(\Sigma_{n+1})}^2)^{1/2}. $$

Further take functions $h, g_1, g_2$ defined on $\mathcal{S}$ as

$$ h = \alpha \| \mathcal{F}_1(\xi, \phi_1, \phi_2, \dots, \phi_{n+1}) - f \|_{\mathcal{D}_1}^2 + \beta \| \mathcal{F}_2(\xi, \phi_1, \phi_2, \dots, \phi_{n+1}) \|_{\mathcal{D}_2}^2, $$

$$ g_1 = \lambda_0 R (\|\xi\|_{[\Xi]^3}), $$

$$ g_2 = \sum_{i=1}^{n+1} \lambda_i \|\phi_i\|_{L^2(\Sigma_i)}^2. $$

Note that the above problem (23) can equivalently be written as the following problem

$$ \min \left\{ H(x_1, x_2) \equiv h(x_1, x_2) + g_1(x_1) + g_2(x_2) \mid (x_1, x_2) \in \mathfrak{B}_1 \times \mathfrak{B}_2 \right\}, $$

where $\mathfrak{B}_1, \mathfrak{B}_2, h, g_1, g_2$ satisfy the following conditions from [32]: (P1)-(P5), (A1) which can easily be shown and

(A2) The partial (Fréchet) derivative of $h$ with respect to the $i$-th component, denoted by $\nabla_i h \in \mathfrak{B}_i^*$, is Lipschitz continuous with Lipschitz constant $L_i \in (0, \infty], i = 1, 2$. with $\min\{L_1, L_2\} < \infty$; exemplarily, for $i = 1$ (analogously for $i = 2$) it holds that $\|\nabla_1 h(x_1 + y_1, x_2) - \nabla_1 h(x_1, x_2)\|_{1,*} \le L_1 \|y_1\|_1$ for all $(x_1, x_2) \in \mathcal{D}, y_1 \in \mathfrak{B}_1$, such that $x_1 + y_1 \in \text{dom } g_1$,

as a consequence of the chain rule of Fréchet differentiation and the fact that $L^p, 1 < p < \infty$ are uniformly convex hence uniformly smooth, see for example [33, Part 3, Chap. II, Sec. 1, Prop. 8 and Part 3, Chap. II, Sec. 2, Prop. 2], hence their norms are uniformly Fréchet differentiable [33, Part 3, Chap. II, Sec. 2, Prop. 1] and that $\mathcal{F}_1, \mathcal{F}_2$ are bounded linear operators hence they are their Fréchet derivatives with respect to $\xi$ and $(\phi_1, \phi_2, \dots, \phi_{n+1})$.
(A3a) Recall that a function $h$ is strongly convex if there exists $\sigma > 0$ such that

$$ \langle \nabla h(x) - \nabla h(\bar{x}), \bar{x} - x \rangle \geq \sigma \|x - \bar{x}\|^2, $$

for all $x \in \text{dom } h$. The function $h: \mathfrak{B}_1 \times \mathfrak{B}_2 \to \mathbb{R}$ is quasi-strongly convex with respect to $\mathcal{O}^*$, with modulus $\sigma > 0$, that is, for all $x \in \mathcal{D}$ and $\bar{x} := \arg \min\{\|x - y\|\|y \in \mathcal{O}^*\}$, the projection of $x$ onto $\mathcal{O}^*$, it holds

$$ h(\bar{x}) \geq h(x) + \langle \nabla h(x), \bar{x} - x \rangle + \frac{\sigma}{2} \|x - \bar{x}\|^2, $$

as a consequence of $h$ being the sum of squares of Hilbert norms each of which is strongly convex function with modulus $\sigma > 0$ and the boundedness of $\mathcal{F}_1, \mathcal{F}_2$ in $\xi$ and $(\phi_1, \phi_2, \dots, \phi_{n+1})$.

We have shown that Problem 1 satisfies the assumption of [32, Thm 1], hence the sequence generated by the alternating minimisation algorithmconverges to the minimum of the Tikhonov functional, $\mathcal{T}_{f,\alpha,\beta,\lambda_0,\dots,\lambda_{n+1}}$, linearly. $\square$

**Remark 4.** For the discrete version of Problem 1 we have that the alternating minimisation sequence converges to a solution of the problem as a consequence of the results provided in [34].

## 5.4 EEG, sEEG and cortical mapping

Numerically, we solve discretised versions of the continuous problems of the previous subsection. Note that in what follows $\xi$ is discretised. In the discrete problem we introduce regularisation matrices $\Gamma_j$ that discretise the regularisers in the continuous problem. For EEG or sEEG we solve the following problem:

**Problem 2.** Given point-wise recordings of electric potential, $\{\phi(Y)\}$, $\alpha, \beta, \lambda_j > 0$, find $(\xi^*, \Phi^*)$ such that

$$ \begin{aligned} (\xi^*, \Phi^*) = \operatorname*{arg\,min}_{\xi \in [\Xi]^3, \Phi \in \mathbb{R}^{m^*}} & \left( \alpha \|\mathbf{F}_E(Y)(\xi, \Phi) - \sigma(Y)\phi(Y)\|_2^2 + \beta \|\mathbf{C}(\xi, \Phi)\|_2^2 \right. \\ & \left. + \lambda_0 R(\|\mathbf{\Gamma}_0 \xi\|_{[\Xi]^3}) + \sum_{j=1}^{n+1} \lambda_j \|\mathbf{\Gamma}_j \Phi\|_2^2 \right), \end{aligned} \quad (25) $$

where $\mathbf{F}_E$ is an appropriate discretisation of the forward model of the electric potential (10) and $\phi(Y)$ is a vector of the potentials at the points $Y = \{y_j\}$ and $m^*$ is the total number of vertices in the triangulation of the cortex, skull and scalp, $\Sigma_{1,T}, \Sigma_{2,T}, \Sigma_{3,T}$, respectively.

## 5.5 MEG and cortical mapping

Numerically, for MEG we will solve the following problem:
**Problem 3.** Given point-wise measurements of components of the magnetic flux density, $\mathbf{B}(Y) \cdot \mathbf{v}(Y)$, $\alpha, \beta, \lambda_j > 0$ find $(\xi^*, \Phi^*)$ such that

$$ (\xi^*, \Phi^*) = \arg\min_{\xi \in [\Xi]^3, \Phi \in \mathbb{R}^{m^*}} \left( \alpha \|\mathbf{F}_B(Y)(\xi, \Phi) \cdot \mathbf{v}(Y) - \mathbf{B}(Y) \cdot \mathbf{v}(Y)\|_2^2 + \beta \|\mathbf{C}(\xi, \Phi)\|_2^2 \right. \\ \left. + \lambda_0 R(\|\mathbf{\Gamma}_0 \xi\|_{[\Xi]^3}) + \sum_{j=1}^{n+1} \lambda_j \|\mathbf{\Gamma}_j \Phi\|_2^2 \right), \quad (26) $$

where $\mathbf{F}_B$ is an appropriate discretisation of the forward model of the magnetic flux density (13) and $\mathbf{B}(Y) \cdot \mathbf{v}(Y)$ is a vector of the $\mathbf{v}$ components of the magnetic flux density at the points $Y = \{y_j\}$ and $m^*$ is the total number of vertices in the triangulation of the cortex, skull and scalp, $\Sigma_{1,T}, \Sigma_{2,T}, \Sigma_{3,T}$, respectively.

## 5.6 EEG, sEEG, MEG and cortical mapping

There is an obvious connection among the EEG, sEEG and MEG source localisation problem and cortical mapping as can be seen in (15 and (20) hence these problems can be solved in a unified way by making appropriates changes to either (25) and (26) if simultaneous recordings of EEG, sEEG and MEG are available.

## 5.7 Resolving the discrete problems

In [34] weaker assumptions than those required in Theorem 2 are stated, in particular the omission of the Fréchet subdifferentiability of $g_1$ and $g_2$ in (P2). This weaker set of assumptions ensures sub-linear convergence of the objective function to the minimum when the Banach spaces $\mathfrak{B}_i$ are $\mathbb{R}^{n_i}$ for some $n_i \in \mathbb{N}$, $i = 1, 2$. This is particularly useful once we have discretised the problem such as in the case measure being discretised, as a collection of dipoles. This discretisation of measures results in the total variation norm of the measures being the $\| \cdot \|_{[t_1(\mathbb{N})]^3}$ of the sequence of Euclidean norms of the dipole moments and this norm is not Fréchet subdifferentiable. This leads to the conclusion that applying an alternating minimisation procedure to the discrete problems for sEEG, EEG and MEG results in obtaining the minimisers that we seek. Note that the solutions we obtain for these problems depend on the regularisation parameters $\alpha, \beta, \lambda_j$. The parameters, $\alpha, \beta, \lambda_j$ are chosen such that the influence of the noise in the measured data on the solution is minimal. Hence, the choice of the parameters $\alpha, \beta, \lambda_j$ is of great importance. Since in practise there may be no information on the properties of the noise corrupting data we propose the use of the L-hypersurface technique for choosing appropriate regularisation parameters, see for example [35], for an exploratory study of this technique. The L-hypersurface approach is a generalisation of the L-curve technique, see for example [36]. In short, in the L-hypersurface technique one aims to find the point of maximum Gaussian curvature on a hypersurface that is generated by plotting the data discrepancy of a solution against the regularisation parameters generating the solution; the plotted values are scaled appropriately. It has to be noted that this a computationally expensive technique and other less computation-
ally expensive techniques can be implemented such as one called the minimal distance function technique which is studied in [35].

# 6 Numerical Results

We now present some numerical results of the inverse source localisation problem using EEG, sEEG and MEG data, we also used combined sEEG and MEG data. These numerical results were obtained using code written in MATLAB and the graphics were produced using a MATLAB add-on Toolbox Graph [37]. The meshes used were processed using the MATLAB add-on Iso2Mesh [38]. Iso2Mesh was used for decimation of the meshes and to fix some defects. The forward problem can be solved satisfactorily with disregard to such defects but the inverse problem is very sensitive to such defects. The processing done in Iso2Mesh resulted in meshes that differed from the ones used in the forward problem thus helping in avoiding the inverse crime. In the forward model the following meshes were used

<table><thead><tr><td></td><td>Number of points</td><td>Number of triangles</td></tr></thead><tbody><tr><td>Grey-white matter interface</td><td>9 312</td><td>18 620</td></tr><tr><td>Cortex</td><td>9 312</td><td>18 620</td></tr><tr><td>Skull</td><td>3 045</td><td>6 078</td></tr><tr><td>Scalp</td><td>3 753</td><td>7 502</td></tr></tbody></table>

The synthetic data we use was generated using OpenMEEG, see [5] and [39], which is based on the symmetric boundary elements method. OpenMEEG uses current dipoles as the elementary electromagnetic object. The dipoles used in OpenMEEG to generate the data were outwardly normally oriented to the grey/white matter interface. In the source recovery problem we attempt to recover the locations of the current dipoles associated with these data by using either $W^{1/2,2}$ functions or a finite collection of dipoles. We attempt the recovery with $[W^{1/2,2}]^3$ vector-field that are normally oriented to the grey/white matter interface or with a collection of dipoles on the barycentres of the triangles of the meshes hence we need only recover the moments of the dipoles placed on the barycentres all triangles. With these source recovery problems we also solve the cortical mapping problem. Note that for the $W^{1/2,2}$ recovery we can use direct inversion methods to solve the inverse problems however for completeness we used the alternating minimisation procedure to solve these problems as well. For solving the inverse problem the following meshes were used for the $[W^{1/2,2}]^3$ source,

<table><thead><tr><td></td><td>Number of points</td><td>Number of triangles</td></tr></thead><tbody><tr><td>Grey-white matter interface</td><td>23 277</td><td>46 550</td></tr><tr><td>Cortex</td><td>27 188</td><td>54 372</td></tr><tr><td>Skull</td><td>3 037</td><td>6 070</td></tr><tr><td>Scalp</td><td>2 416</td><td>4 828</td></tr></tbody></table>

where as for the dipolar source the grey-white matter interface had 13 956 triangles and 6 980 points. The recovery for the synthetic data was done with no noise added. We performed 1 000 iterations of the alternating minimisation procedure. For the dipolar
case the FISTA step was done with a maximum of 1 000 iterations. The computation
time was in the order of hours for the [$W^{1/2,2}]^3$ recovery and for the dipolar recovery
the computation time was more than 24 hours but less than 48 hours.

Figures 3, 4, 5, 6, 7, 8, 9 and 10 show the ground truth (OpenMEEG data) and recovery results obtained from using sEEG, MEG, combined sEEG and MEG and EEG data, respectively. In each figure the first column shows the ground truth data generated by OpenMEEG using the dipoles represented as red dots in the bottom most figure of the first column. The second column of each Figure shows the recovery results when the source is assumed to be [$W^{1/2,2}]^3$ vector-field normally oriented to the grey/white matter interface and the third column shows the recovery results when the source is assumed to a finite collection of dipoles, the dots representing the dipoles are colour coded by the magnitude of their moments. In the case of [$W^{1/2,2}]^3$ sources, the source is interpreted as being located within a neighbourhood of the various local maxima of the recovered $W^{1/2,2}$ function supported on the grey/white matter interface. For dipoles, the area with the highest concentration of dipoles with moments of largest magnitude is taken to be the area supporting the true source.

From Figures 3 and 7 we see that the placement of sEEG electrodes has a higher
influence on the source recovery compared to EEG and MEG. The figures suggest
that the closer the sEEG electrodes are to the true source the better the recovery
of the source. We also observe that the recovered surface electrical potentials are
overestimated on the skull and the scalp while being underestimated on the cortex
however their distribution is in accordance with the recovered source. It seems that
the sharper the recovered source the more accurate the recovered surface electrical
potentials in terms of their distribution of the surfaces.

In Figures 4 and 8 we immediately note that the surface electrical potential are underestimated however their distributions are more representative of the ground truth compared to the recovery managed using sEEG data. With MEG data the localisation of the sources is closer to the true sources and also more spatially localised than the ones recovered with sEEG data when the sEEG electrodes are far from the true source. This suggests that for source localisation MEG is more robust than sEEG. That we managed to recover surfaces electric potentials that have the correct distribution with MEG data is short of impressive, especially when look at the recovery of the surface electric potentials with a [$W^{1/2,2}]^3$ source in both figures.

Observing that a better source localisation leads to be a better recovery of the surface
electric potentials, it is of interest to see if the robustness of the source localisation of
MEG can be complimented with the electric data of sEEG for a better recovery of the
source and the surface electric potentials. Figures 5 and 9 show results obtained by
using simultaneous MEG and sEEG data. We can observe that the recovered source
and the recovered surface electrical potentials have properties that are shared between
the solutions from each modality hence represent a better and more robust recovery.
The source recovered with combined MEG and sEEG data is more spatially localised
than with standalone modalities and it is located closer to the true source. Further,
the recovered surface potential are much more representative of the ground truth than
had been previously seen with standalone modalities.

It is evident from Figures 6 and 10 that the recovery done with EEG data outper-
forms the other recoveries. The method we implemented here makes the recovery with EEG data robust as it managed a near perfect recovery of both sources and surface electrical potential. In Figure 10 we see that the source recovered for [$W^{1/2,2}]^3$ is not as spatially localised as the dipolar source, however, with the interpretation of source we are using, the source location of the [$W^{1/2,2}]^3$ overlaps with the true source location. We think that the “depth” of the source in this case resulted in the poor spatial location we observed here.

We applied the alternating minimisation algorithm to real auditory data on the same meshes we used for the synthetic data. The results are shown in Figures 11 and 12. The data is from the supplementary material of the paper [3] which had been cleaned and filtered using independent component analysis using EEGLAB. In our recovery of the source we interpret the result as to show a main activation of the region around corpus callosum in blue. When looking at the left and right sides skull and scalp there are local peaks in the potential in dark red and we observe the same for the electric potential on the left and right lateral sulci of the cortex. We thus conclude that there are secondary activations in the “sulci” of the grey white matter that correspond to the lateral sulci of the cortex. These secondary activations agree with recovery results presented in [3] for the N100 auditory evoked response. The manner in which we made the analysis here demonstrates some of the benefits of performing the cortical mapping and source localisation simultaneous as we managed to trace the secondary activation from peaks of the electric potential on the scalp. This additional information helps to reinforce the conclusions of the source recovery.

What is evident from the numerical results presented here is that the surface electrical potential on the cortex was poorly recovered across all the modalities. We put this to the complicated geometry of the surface in which it folds on itself. Even though we could not show it in the graphics, the maxima of the surface electric potential on the cortex are located on the “sulci” of the cortex in these examples. We managed to recreate these maxima in approximately the correct locations on the sulci of the cortex.
<figure><img src="image_4.png" /><figcaption>Figure 2: Above is the head model used in the computations, the model had 4 surfaces, the inner most red surface represents the grey/white matter interface, the dark blue surface represents the cortical surface, the light blue surface represents the outer surface of the skull and the outer most yellow surface represents the scalp. In the same figure we show the location of the 198 sEEG electrodes as blue dots, the majority of them being in the region enclosed by the grey/white matter interface. The 64 EEG electrodes are represented as green dots on the scalp and the locations of the 151 Superconducting Quantum Interference Devices (SQUIDs) for recording MEG measurements are shown as black dots outside the head.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Recovery with sEEG data.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Recovery with MEG data.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Recovery with combined MEG and sEEG data.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Recovery with EEG data.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Recovery with SEEG data.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Recovery with MEG data.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Recovery with combined MEG and sEEG data.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: Recovery with EEG data.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11: [$W^{1/2,2}]^3$ recovery with real auditory EEG data for the N100 response from supplementary data of the paper [3].</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12: Dipole source recovery with real auditory EEG data for the N100 response from supplementary material of the paper [3].</figcaption></figure>
# 7 Conclusions

We have presented a method that uses the double and single layer potential operators only in the transmission of the electric potential and the magnetic flux density associated with a source which in turn allows to solve the inverse source recovery problem and the inverse cortical mapping problem using either electric potential data or magnetic flux density data. The simultaneous resolution of these problems allow to capture the full behaviour of the source which possibly aids in having better solutions for both inverse problems. The method is also less computationally complex as it involves fewer boundary integral operators than the symmetric boundary elements method that has been used previously in inverse cortical mapping. The method also allows the use of realistic geometries of the head and making exact computations at arbitrary points in space which allows for exact placement of sensors in the models by taking advantage of the formulas of [23] for improved numerical accuracy for building and solving the discretised problems. The manner in which the discretised problems are built allows for the use of numerical frames/bases of Banach spaces defined on surfaces hence allowing the numerical implementation of the method to the Banach spaces in which the problem is solvable as was demonstrated in Theorem 1. This opens up the possibility for models that are closer to reality than have been previously done. The alternating minimisation procedure employed allows the use of discretisations that offer the best accuracy for the source coupled with the freedom to use the best and efficient methods for recovering the source and associated surface potentials. Further the method offers a natural coupling of the electric data and magnetic flux density making it easier to combine these data for the recovery as demonstrated in the combined use of sEEG and MEG data.

Future work on this subject includes making performance comparisons with existing methods for source recovery and inverse cortical mapping, testing the method on real data and implementing this method on a wider range of source classes that encourage the recovery of certain assumed behaviours of brain activities.

# 8 Acknowledgements

The work of M. Nemaire is supported by the Agence Nationale de la Recherche (ANR) through the the project ANR-18-CE40-0035.

The authors are grateful to the OPAL infrastructure from Université Côte d'Azur for providing resources and support.

<footnote n="1">¹ FACTAS, Inria Sophia Antipolis-Méditerranée, 2004 Rte des Lucioles, 06902 Valbonne, France</footnote>

<footnote n="2">² Institut de Neurosciences des Systèmes, Aix-Marseille Université, INS/Inserm, 13005 Marseille, France</footnote>

<footnote n="3">³ Center of Applied Mathematics, Ecole des Mines ParisTech, CS 10 207, 06904 Sophia Antipolis Cedex, France</footnote>
⁴ Institut de Mathématiques de Bordeaux UMR 5251, Université de Bordeaux, 351, cours de la Libération, 33405 Talence, France

## References

[1] A. Gramfort, M. Luessi, E. Larson, D. Engemann, D. Strohmeier, C. Brodbeck, L. Parkkonen, and M. Hämäläinen. MNE software for processing MEG and EEG data. *NeuroImage*, 86, 10 2013.

[2] M. Clerc, J. Leblond, J.-P. Marmorat, and T. Papadopoulo. Source localization using rational approximation on plane sections. *Inverse Problems*, 28, 05 2012.

[3] M. Stropahl, A.-K. Bauer, S. Debener, and M. Bleichner. Source-Modeling Audiory Processes of EEG Data Using EEGLAB and Brainstorm. *Frontiers in Neuroscience*, 12, 05 2018.

[4] M. Clerc and J. Kybic. Cortical mapping by Laplace-Cauchy transmission using a boundary element method. *Inverse Problems*, 23:2589-2601, 12 2007.

[5] J. Kybic, M. Clerc, T. Abboud, O. Faugeras, R. Keriven, and T. Papadopoulo. A common formalism for the integral formulations of the forward EEG problem. *IEEE transactions on medical imaging*, 24:12-28, 02 2005.

[6] L. Rahmouni, S. Adrian, K. Cools, and F. Andriulli. Conforming Discretizations of Boundary Element Solutions of the Electroencephalography Forward Problem. *Comptes Rendus Physique*, 19:7-25, 01 2018.

[7] V.D. Kupradze and M.A. Aleksidze. The method of functional equations for the approximate solution of certain boundary value problems. *USSR Computational Mathematics and Mathematical Physics*, 4(4):82-126, 1964.

[8] Liang Yan, Chu-Li Fu, and Feng-Lian Yang. The method of fundamental solutions for the inverse heat source problem. *Engineering Analysis with Boundary Elements*, 32(3):216-222, 2008.

[9] D. Cosandier-Rimélé, J.-M. Badier, P. Chauvel, and F. Wendling. A Physiologically Plausible Spatio-Temporal Model for EEG Signals Recorded With Intracerebral Electrodes in Human Partial Epilepsy. *IEEE transactions on bio-medical engineering*, 54:380-8, 04 2007.

[10] J. Sarvas. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. *Physics in medicine and biology*, 32 1:11-22, 1987.

[11] D. Shattuck and R. Leahy. Brainsuite: An Automated Cortical Surface Identification Tool. *Medical Image Analysis*, 6:129-42, 07 2002.
[12] R. Ilmoniemi and J. Sarvas. *Brain Signals: Physics and Mathematics of MEG and EEG*. The MIT Press, Cambridge, MA, 01 2019.

[13] M. Hämäläinen, R. Hari, R. Ilmoniemi, J. Knuutila, and O. Lounasmaa. Magnetoencephalography: Theory, instrumentation, and applications to noninvasive studies of the working human brain. *Rev. Mod. Phys.*, 65:413–, 04 1993.

[14] P. de Lange, E. Boto, N. Holmes, R. M. Hill, R. Bowtell, V. Wens, X. De Tiège, M. J. Brookes, and M. Bourguignon. Measuring the cortical tracking of speech with optically-pumped magnetometers. *NeuroImage*, 233:117969, 2021.

[15] G. Verchota. Layer potentials and regularity for the Dirichlet problem for Laplace's equation in Lipschitz domains. *Journal of Functional Analysis*, 59(3):572–611, 1984.

[16] R. Kress. *Linear Integral Equations*. Applied Mathematical Sciences. Springer-Verlag Berlin Heidelberg, 1989.

[17] Eugene Fabes, Osvaldo Mendez, and Marius Mitrea. Boundary Layers on Sobolev-Besov Spaces and Poisson's Equation for the Laplacian in Lipschitz Domains. *Journal of Functional Analysis*, 159(2):323–368, 1998.

[18] S. Axler, P. Bourdon, and W. Ramey. *Harmonic Function Theory*. Graduate Texts in Mathematics. Springer New York, 2006.

[19] J. D. Jackson. *Classical electrodynamics*. Wiley, New York, NY, 1962.

[20] R. Dautray and Lions J.-L. *Mathematical Analysis and Numerical Methods for Science and Technology, Volume 1: Physical Origins and Classical Methods*. Springer-Verlag Berlin Heidelberg, 2000.

[21] E. Corona and S. Veerapaneni. Boundary integral equation analysis for suspension of spheres in Stokes flow. *Journal of Computational Physics*, 362:327–345, 2018.

[22] J. Rahol and S. Tissari. Iterative solution of dense linear systems arising from the electrostatic integral equation in MEG. *Physics in medicine and biology*, 47:961–75, 04 2002.

[23] S. Nintcheu Fata. Explicit expressions for 3D boundary integrals in potential theory. *International Journal for Numerical Methods in Engineering*, 78:32–47, 2009.

[24] D. Mitrea, M. Mitrea, and J. Pipher. Vector potential theory on nonsmooth domains in $R^3$ and applications to electromagnetic scattering. *Journal of Fourier Analysis and Applications*, 3:131–192, 1997.

[25] L. Baratchart, C. Villalobos Guillen, D. P. Hardin, M. C. Northington, and E. B. Saff. Inverse Potential Problems for Divergence of Measures with Total Variation Regularization. *Foundations of Computational Mathematics*, 20:1273–1307, 2020.
[26] S. K. Smirnov. Decomposition of solenoidal vector charges into elementary solenoids, and the structure of normal one-dimensional flows. *St. Petersburg Math. J.*, 5:841–867, 1994.

[27] R. A. Adams. *Sobolev spaces* / Robert A. Adams. Academic Press New York, 1975.

[28] F. Demengel, G. Demengel, and R. Erné. *Functional Spaces for the Theory of Elliptic Partial Differential Equations*. Universitext. Springer London, 2012.

[29] B. Hofmann, B. Kaltenbacher, C. Pöschl, and O Scherzer. A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators. *Inverse Problems*, 23:987, 04 2007.

[30] H. Brezis. *Functional Analysis, Sobolev Spaces and Partial Differential Equations*. Universitext. Springer New York, 2010.

[31] J. Bolte, S. Sabach, and M. Teboulle. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. *Mathematical Programming*, 146:459–494, 08 2013.

[32] J. Both. On the rate of convergence of alternating minimization for non-smooth non-strongly convex optimization in Banach spaces. *Optimization Letters*, pages 1–15, 05 2021.

[33] B. Beauzamy. *Introduction to Banach spaces and their geometry*. Number no. 86; no. 1985 in Notas de matemática. North-Holland, 1985.

[34] A. Beck. On the Convergence of Alternating Minimization for Convex Programming with Applications to Iteratively Reweighted Least Squares and Decomposition Schemes. *SIAM J. Optim.*, 25:185–209, 2015.

[35] M. Belge, M. Kilmer, and E. Miller. Efficient Determination Of Multiple Regularization Parameters In A Generalized L-Curve Framework. *Inverse Problems*, 18:1161–1183, 09 1999.

[36] P. C. Hansen. The L-Curve and its Use in the Numerical Treatment of Inverse Problems. In *Computational Inverse Problems in Electrocardiology*, ed. P. Johnston, *Advances in Computational Bioengineering*, pages 119–142. WIT Press, 2000.

[37] G. Peyre. Toolbox Graph, 2020.

[38] Q. Fang. Iso2mesh, 2020.

[39] A. Gramfort, T. Papadopoulo, E. Olivi, and M. Clerc. OpenMEEG: opensource software for quasistatic bioelectromagnetics. *BioMedical Engineering OnLine*, 9, 2010.