Algebraic Reconstruction of Current Dipoles and Quadrupoles in Three-Dimensional Space

Abstract

This paper presents an algebraic method for an inverse source problem for the Poisson equation where the source consists of dipoles and quadrupoles. This source model is significant in the magnetoencephalography inverse problem. The proposed method identifies the source parameters directly and algebraically using data without requiring an initial parameter estimate or iterative computation of the forward solution. The obtained parameters could be used for the initial solution in an optimization-based algorithm for further refinement.

1. Introduction

Inverse source problems for the Poisson equation have a wide variety of applications such as bioelectromagnetic inverse problems. Magnetoencephalography (MEG) and electroencephalography (EEG) are typical examples in which the current source inside the human brain is inversely reconstructed from measurements of the magnetic field outside the head and the electric potential on the head surface [1–3]. Due to the quasistatic property of the magnetic and electric fields, MEG and EEG exhibit high temporal resolution compared to other tools for the noninvasive inspection of brains. Hence, the development of a stable inverse algorithm which reconstructs the source with high spatial resolution is of importance for MEG and EEG.

Conventional algorithms for the MEG inverse problem, whose solution is generally nonunique [4], can be divided into two categories [2]: the parametric and imaging approaches. The parametric approach assumes that the current source inside the brain is expressed by a relatively small number of equivalent current dipoles (ECDs), which guarantees the uniqueness of the solution [4] and finds their number, positions, and dipole moments. A typical method uses optimization-based algorithms to minimize the squared sum of the difference between the forward solution and data measured at a finite number of sensor positions (see, e.g., [5]). The forward solution is given by either the Geselowitz formula, which assumes that the head consists of arbitrarliy-shaped layered domains, or the Sarvas formula, which assumes that the head consists of concentric spheres [1]. A problem with optimization-based algorithms is that they require an initial solution close to the true solution, or the algorithm will often converge to a local minimum. To resolve this problem, algebraic methods have been proposed and are gaining increased importance [6–10]. In these methods, the ECD parameters are reconstructed directly and algebraically from the boundary measurements of the magnetic field.

In contrast, the second category of inverse algorithms, the imaging approach, assumes that the many dipoles are set at the vertices of a mesh overlaying the cortical surface and solves for their dipole moments. An advantage of this approach is that the spatially distributed source can be reconstructed by solving linear equations for the unknown dipole moments. However, it has the disadvantages that the solution is not unique and an oversmoothed solution is often obtained by regularization, such as the minimum 𝐿₂ norm solution [11]. A method using the minimum 𝐿₁ norm solution has also been proposed to reconstruct a sparse source [12].

Recently, a new method, the multipolar representation of the source, has been developed that incorporates some of the advantages of the above two methods, and has attracted considerable attention [13–18]. In this model, instead of expressing the current source by an equivalent current dipole, an equivalent dipole and quadrupole [15, 16] or an equivalent dipole and octopole [13, 14] are used where the quadrupole or octopole is determined depending on the spatial extent of the support of the current source. Jerbi et al. showed that the centroid of the spatially distributed source, which they called a patch source, can be estimated more accurately using the dipole and quadrupole model than using the dipole model by means of the nonlinear least squares method [16]. We considered a two-dimensional (2D) problem using a complex analysis framework and proposed an algebraic method to reconstruct the dipole and quadrupole parameters directly from the boundary measurements of the meromorphic function [17]. The aim of the current paper is to extend our algebraic method from the 2D case to the 3D case so that the dipole and quadrupoles, which equivalently represent the neural current, can be reconstructed from the magnetic field data without an initial parameter estimate or iterative computing of the forward solution. In [19, 20], we proposed an algebraic method when the dipoles were distributed in a plane parallel to the 𝑥𝑦-plane, which is a very special case and is severely restricted in its practical usage. This paper derives a method for the general case. From a practical viewpoint, our method can provide an estimate of the number of patches as well as a good initial solution close to the true solution for optimization-based algorithms.

This paper is organized as follows. In Section 2, the forward problem with the dipole-quadrupole source model is summarized and the inverse problem is formulated. In Section 3, we propose a method to reconstruct the 𝑥𝑦-coordinates of the dipole-quadrupole source by solving simultaneous algebraic equations of second degree. A method to estimate the 𝑧-coordinates is also proposed. Section 4 is devoted to numerical simulations for estimating spatially distributed dipoles.

2. Forward and Inverse Problems with the Dipole and Quadrupole Source Model

Let Ω₁, Ω₂, Ω₃, and Ω₄ be concentric balls centered at the origin in 3D space, where Ω_𝑖⊂Ω_𝑖+1 for 𝑖=1,2,3, as shown in Figure 1. Here, Ω₁, Ω₂/Ω₁, and Ω₃/Ω₂ represent the brain, skull, and scalp, respectively. Ω₃ represents the head. We assume that the radial component of the magnetic field is measured on the sphere Γ=𝜕Ω₄ with the radius of 𝑅. Although we use this simple head model as well as the spherical sensor surface in this paper, our method can be extended to a more realistic case when the head is modeled by a piecewise homogeneous layered domain and the sensors are set on an arbitrarily shaped surface as assumed in the algebraic method for the dipole source model [10]. Let 𝐉_𝑝(𝐫) represent the neural currents in the innermost domain Ω₁. The magnetic permeability is assumed to be constant 𝜇₀ in the whole space.

Figure 1 

Head and sensor model.

First, we derive the representation of the radial magnetic field in terms of the equivalent dipoles and quadrupoles, which has already been given in [15], to determine how the quadrupole terms equivalently express the spatial extent of the source.

It is known that with this head model, the radial component of the magnetic field at 𝐫 is equal to that when 𝐉_𝑝(𝐫) exists in infinite homogeneous space and is given by the Biot-Savart law: (𝐫)=𝜇𝐫⋅𝐁⁰4𝜋_Ω1∇^′1𝐫−𝐫^′×𝐫⋅𝐉_𝑝𝐫^′𝑑𝑣^′. (1) Let us assume that the neural activity is localized at several domains 𝐷_𝑘(⊂Ω₁) for 𝑘=1,2,…,𝑁 so that 𝐉_𝑝(𝐫) is given by 𝐉_𝑝(𝐫)=^𝑁_𝑘=1𝐣_𝑝(𝐫)𝜒_𝐷𝑘, (2) where 𝜒_𝐷𝑘 is the characteristic function. Then, the radial component of the magnetic field is given by (𝐫)=𝜇𝐫⋅𝐁⁰4𝜋^𝑁_{𝑘=1𝐷𝑘}∇^′1𝐫−𝐫^′×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′. (3) Assume that 𝐷_𝑘 is contained in a ball 𝐵_𝑘⊂Ω₁, and there exists a “centroid” 𝐫_𝑘 of 𝐷_𝑘 such that 𝐫−𝐫_𝑘>𝐫^′−𝐫_𝑘,∀𝐫∈Γ,𝐫^′∈𝐷_𝑘. (4) From the Taylor series expansion, we have ∇^′1𝐫−𝐫^′=∇^′1𝐫−𝐫^′_{𝐫′=𝐫𝑘}+∇^′∇^′1𝐫−𝐫^′_{𝐫′=𝐫𝑘}𝐫^′−𝐫_𝑘+⋯. (5) Substituting this expansion into the right-hand side of (3) gives (𝐫)=𝜇𝐫⋅𝐁⁰4𝜋^𝑁_𝑘=1𝐫−𝐫_𝑘𝐫−𝐫_𝑘³⋅𝐫_𝑘×_𝐷𝑘𝐣_𝑝𝐫^′𝑑𝑣^′+𝜇₀4𝜋^𝑁_𝑘=13𝐫×𝐫_𝑘𝐫−𝐫_𝑘+𝐫−𝐫_𝑘²𝑋_𝑟𝐫−𝐫_𝑘⁵∶_𝐷𝑘𝐣_𝑝𝐫^′𝐫^′−𝐫_𝑘𝑑𝑣^′+⋯, (6) where 𝑋_𝑟 is the cross-product tensor defined by 𝑋_𝑟𝐚=𝐫×𝐚, (7) and hence is written as 𝑋_𝑟=0−𝑧𝑦𝑧0−𝑥−𝑦𝑥0for𝐫=𝑥,𝑦,𝑧^𝑇. (8) In (6), “∶” represents the tensor contraction defined by 𝐴∶𝐵=∑³_𝑖=1∑³_𝑗=1𝐴_𝑖𝑗𝐵_𝑖𝑗 for second-order tensors 𝐴 and 𝐵 whose (𝑖,𝑗) components are given by 𝐴_𝑖𝑗 and 𝐵_𝑖𝑗, respectively. See Appendix A for the derivation of (6).

Now, let 𝐩_𝑘≡_𝐷𝑘𝐣_𝑝𝐫^′𝑑𝑣^′,𝑄_𝑘≡12_𝐷𝑘𝐣_𝑝𝐫^′𝐫^′−𝐫_𝑘+𝐫^′−𝐫_𝑘𝐣_𝑝𝐫^′𝑑𝑣^′=𝑞_{𝑥𝑥,𝑘}𝑞_{𝑥𝑦,𝑘}𝑞_{𝑥𝑧,𝑘}𝑞_{𝑦𝑥,𝑘}𝑞_{𝑦𝑦,𝑘}𝑞_{𝑦𝑧,𝑘}𝑞_{𝑧𝑥,𝑘}𝑞_{𝑧𝑦,𝑘}𝑞_{𝑧𝑧,𝑘}, (9) then we have (𝐫)=𝜇𝐫⋅𝐁⁰4𝜋^𝑁_𝑘=1𝐫_𝑘×𝐩_𝑘⋅𝐫−𝐫_𝑘𝐫−𝐫_𝑘³+𝜇₀4𝜋^𝑁_𝑘=1𝑄_𝑘∶3𝐫×𝐫_𝑘𝐫−𝐫_𝑘+𝐫−𝐫_𝑘²𝑋_𝑟𝐫−𝐫_𝑘⁵+⋯. (10) Note here that we define 𝑄_𝑘 to be symmetric, 𝑞_{𝑖𝑗,𝑘}=𝑞_{𝑗𝑖,𝑘}, since in (6) only the symmetric part of ∫_𝐷𝑘𝐣_𝑝(𝐫^′)(𝐫^′−𝐫_𝑘)𝑑𝑣^′ contributes to 𝐫⋅𝐁 as shown in [15]. Equation (10) was given in (46) in [15] (when 𝑁=1). The first term in (10) is the magnetic field created by an “equivalent current dipole” 𝐩_𝑘, and the second term is the magnetic field created by an “equivalent current quadrupole” 𝑄_𝑘. Note here that 𝐩_𝑘 does not depend on the size of 𝐷_𝑘, while 𝑄_𝑘 depends on the spatial extent of 𝐣_𝑝 in 𝐷_𝑘. 𝑄_𝑘 is a 3×3 tensor of order 2 and is called the quadrupole moment tensor.

Here, we can show that the truncation of (10) up to the second term, (𝐫)=𝜇𝐫⋅𝐁⁰4𝜋^𝑁_𝑘=1𝐫_𝑘×𝐩_𝑘⋅𝐫−𝐫_𝑘𝐫−𝐫_𝑘³+𝑄_𝑘∶3𝐫×𝐫_𝑘𝐫−𝐫_𝑘+𝐫−𝐫_𝑘²𝑋_𝑟𝐫−𝐫_𝑘⁵, (11) is the forward solution for the source model given by 𝐉_𝑝(𝐫)=^𝑁_𝑘=1𝐩_𝑘𝛿𝐫−𝐫_𝑘+^𝑁_𝑘=1𝑄_𝑘∇𝛿𝐫−𝐫_𝑘, (12) where 𝛿 is the Dirac delta function (see Appendix B). We call (12) the “equivalent dipole and quadrupole source model” or simply the “dipole-quadrupole source model.” We are now ready to state our inverse problem.

Problem 1. Assume that a head Ω₃ is modeled by a spherically layered domain. Assume that the neural current is expressed by (12) supported in the inner-most ball Ω₁ representing the brain. Then, reconstruct the number 𝑁 and positions 𝐫_𝑘 of the dipole-quadrupole source in (12) from the radial component of the magnetic field measured on the spherical boundary Γ which encloses Ω₃.

Remark 1. Once the number and positions of the dipole-quadrupole source are determined, the dipole moments 𝐩_𝑘 and quadrupole moment tensors 𝑄_𝑘 can be obtained by solving the relevant linear equation (11). In this paper, we concentrate on reconstructing 𝑁 and 𝐫_𝑘.

3. Algebraic Reconstruction Method

According to (88) in [15], we can express (1) by the multipole expansion 𝐫⋅𝐁(𝐫)=𝜇₀^∞𝑙_{𝑙=0𝑚=−𝑙}𝑙+1𝑀2𝑙+1^𝑙𝑚𝑌^∗_𝑙𝑚𝜃,𝜙𝑟^𝑙+1, (13) where 𝑌_𝑙𝑚𝜃,𝜙=2𝑙+1()!4𝜋𝑙−𝑚(𝑙+𝑚)!×𝑌_𝑙𝑚𝑌𝜃,𝜙∶normalizedsphericalharmonics,_𝑙𝑚𝜃,𝜙=𝑃^𝑚_𝑙(cos𝜃)e^𝑖𝑚𝜙∶sphericalharmonics, (14) and 𝑃^𝑚_𝑙(cos𝜃) are the associated Legendre polynomials. As shown in (84) in [15], the multipole moment has the following relationship with 𝐉_𝑝: 𝑀_𝑙𝑚=1𝑙+1_Ω1∇^′𝑟^′𝑙𝑌_𝑙𝑚𝜃^′,𝜙^′⋅𝐫^′×𝐉_𝑝𝐫^′𝑑³𝑟^′. (15) Substituting (12) into (15) and using (B.1) gives 𝑀_𝑙𝑚=1𝑙+1^𝑁_𝑘=1∇^′𝑟^′𝑙𝑌_𝑙𝑚𝜃^′,𝜙^′_{𝐫′=𝐫𝑘}+𝑄_𝑘∶∇^′∇^′𝑟^′𝑙𝑌_𝑙𝑚𝜃^′,𝜙^′×𝐫^′𝑇_{𝐫′=𝐫𝑘}. (16) On the other hand, using the orthogonality of the spherical harmonics, the multipole moment is shown to have another relationship with the radial magnetic field on a spherical surface Γ [21]: 𝑀_𝑙𝑚=2𝑙+1𝜇₀(𝑙+1)_Γ𝐧^′⋅𝐁𝐫^′𝑅^𝑙𝑌_𝑙𝑚𝜃^′,𝜙^′𝑑𝑆^′, (17) where 𝐧^′ is the outward unit normal vector at 𝐫^′ on Γ. Equating (16) and (17) for 𝑙≥𝑚≥0 gives algebraic equations relating the unknown parameters of the sources to the radial magnetic field on Γ. Although (17) requires the radial components of the magnetic field on the whole boundary Γ, Taulu et al. proposed a new approach for obtaining 𝑀_𝑙𝑚 from the data on the part of Γ by using the truncated spherical harmonic expansion [22]. Using the method, we have algebraic equations about the unknown source parameters with the practically measurable MEG data.

3.1. Algebraic Reconstruction of the 𝑥- and 𝑦-Coordinates

As we derived the algebraic method for the dipole source model in [10], we now use the multipole moments of the sectorial harmonics represented by 𝑙=𝑚=𝑛, where 𝑛∈ℕ. From the relationship 𝑟^𝑛𝑌_𝑛𝑛(𝜃,𝜙)=(2𝑛−1)!!(𝑥+𝑖𝑦)^𝑛 [23], we have ^𝑁_𝑘=1∇𝑥+𝑖𝑦^𝑛_{𝐫=𝐫𝑘}⋅𝐫_𝑘×𝐩_𝑘+𝑄_𝑘∶∇∇𝑥+𝑖𝑦^𝑛×𝐫^𝑇_{𝐫=𝐫𝑘}=2𝑛+1𝜇_0Γ𝐧⋅𝐁(𝐫)𝑥+𝑖𝑦^𝑛𝑑𝑆. (18) The primes in both sides are removed for simplicity. Since it holds that∇𝑥+𝑖𝑦^𝑛1𝑖0=𝑛𝑥+𝑖𝑦^𝑛−1,∇∇𝑥+𝑖𝑦^𝑛×𝐫^𝑇𝑖()𝑧=𝑛𝑛−1𝑥+𝑖𝑦^𝑛−2−(𝑛−1)𝑧𝑥+𝑖𝑦^𝑛−2𝑖𝑥+𝑖𝑦^𝑛−1−(𝑛−1)𝑧𝑥+𝑖𝑦^𝑛−2−𝑖(𝑛−1)𝑧𝑥+𝑖𝑦^𝑛−2−𝑥+𝑖𝑦^𝑛−1−𝑖𝑛𝑥+𝑖𝑦^𝑛−1𝑛𝑥+𝑖𝑦^𝑛−10, (19)we obtain ∇𝑥+𝑖𝑦_{𝐫=𝐫𝑘}⋅𝐫_𝑘×𝐩_𝑘=𝑛𝐫_𝑘×𝐩_{𝑘𝑥+𝑖𝑦}𝑥_𝑘+𝑖𝑦_𝑘^𝑛−1,𝑄_𝑘∶∇∇𝑥+𝑖𝑦^𝑛×𝐫^𝑇_{𝐫=𝐫𝑘}=𝑛(𝑛−1)𝑖𝑞_{𝑥𝑥,𝑘}−𝑞_{𝑦𝑦,𝑘}−2𝑞_{𝑥𝑦,𝑘}𝑧_𝑘+𝑞_{𝑦𝑧,𝑘}−𝑖𝑞_{𝑥𝑧,𝑘}𝑥_𝑘+𝑖𝑦_𝑘𝑥_𝑘+𝑖𝑦_𝑘^𝑛−2, (20)(21) where [𝐫_𝑘×𝐩_𝑘]_{𝑥+𝑖𝑦} represents (𝑥-component of 𝐫_𝑘×𝐩_𝑘) + 𝑖×(𝑦-component of 𝐫_𝑘×𝐩_𝑘). See Appendix C for the derivation of (21). Substituting (20) and (21) into (18) gives ^𝑁_𝑘=1𝐫_𝑘×𝐩_{𝑘𝑥+𝑖𝑦}𝑥_𝑘+𝑖𝑦_𝑘^𝑛−1+()×𝑖𝑞𝑛−1_{𝑥𝑥,𝑘}−𝑞_{𝑦𝑦,𝑘}−2𝑞_{𝑥𝑦,𝑘}𝑧_𝑘+𝑞_{𝑦𝑧,𝑘}−𝑖𝑞_{𝑥𝑧,𝑘}𝑥_𝑘+𝑖𝑦_𝑘𝑥_𝑘+𝑖𝑦_𝑘^𝑛−2=2𝑛+1𝜇₀𝑛_Γ𝐧⋅𝐁𝑥+𝑖𝑦^𝑛𝑑𝑆. (22) We now put 𝜁_𝑘=𝑥_𝑘+𝑖𝑦_𝑘,𝜇_𝑘=𝐫_𝑘×𝐩_{𝑘𝑥+𝑖𝑦},𝜈_𝑘=𝑖𝑞_{𝑥𝑥,𝑘}−𝑞_{𝑦𝑦,𝑘}−2𝑞_{𝑥𝑦,𝑘}𝑧_𝑘+𝑞_{𝑦𝑧,𝑘}−𝑖𝑞_{𝑥𝑧,𝑘}𝑥_𝑘+𝑖𝑦_𝑘,𝑐_𝑛=2𝑛+3𝜇₀(𝑛+1)_Γ𝐧⋅𝐁𝑥+𝑖𝑦^𝑛+1𝑑𝑆 (23) and rewrite (22) as ^𝑁_𝑘=1𝜇_𝑘𝜁^𝑛_𝑘+𝑛^𝑁_𝑘=1𝜈_𝑘𝜁_𝑘^𝑛−1=𝑐_𝑛𝑛∈ℕ^∗. (24) Comparing with the special case where the dipoles were distributed in a plane parallel to the 𝑥𝑦-plane in [20], one finds that 𝜈_𝑘 has an extra term: (𝑞_{𝑦𝑧,𝑘}−𝑖𝑞_{𝑥𝑧,𝑘})(𝑥_𝑘+𝑖𝑦_𝑘).

Equation (24) has the same form as that of (6) in [17], and consequently the dipole-quadrupole position projected on the 𝑥𝑦-plane, 𝜁_𝑘(𝑘=1,2,…,𝑁), can be algebraically reconstructed from 𝑐_𝑛(𝑛=0,1,…,2𝑁−1), which is shown as follows. In Section 3.1, we assume that 𝑁 is known a priori and construct an algorithm to estimate 𝜁_𝑘 in (24). A method to estimate 𝑁 as well as 𝜇_𝑘 and 𝜈_𝑘 is then described in Section 3.2.

Let us define 𝜎₁,𝜎₂,…,𝜎_𝑁 by the coefficients of ^𝑁_𝑘=1𝜁−𝜁_𝑘=𝜁^𝑁+𝜎₁𝜁^𝑁−1+⋯+𝜎_𝑁. (25) It was proved in [17] that (24) can be transformed into the second-degree equations for 𝜎₁,𝜎₂,…,𝜎_𝑁: 𝝈^𝑇𝐻_𝑁,𝑚𝝈=0𝑚∈ℕ^∗, (26) where 𝝈=(𝜎_𝑁,…,𝜎₁,1)^𝑇 and 𝐻_𝑁,𝑚 is the (𝑁+1)×(𝑁+1) Hankel matrix given by 𝐻_𝑁,𝑚=𝑐_𝑚𝑐_𝑚+1⋯𝑐_𝑚+𝑁𝑐_𝑚+1𝑐_𝑚+2⋯𝑐_𝑚+𝑁+1𝑐⋮⋮⋮⋮_𝑚+𝑁𝑐_𝑚+𝑁+1⋯𝑐_𝑚+2𝑁. (27) Moreover, it has been shown that the 2𝑁 equations (26) for 𝑚=0,1,…,2𝑁−1 can be turned into linear equations as follows. First, from those 2𝑁 equations, we make 2𝑁−1 equations 𝝈^𝑇𝐻⁽¹⁾_𝑁,𝑚𝝈=0(𝑚=0,1,…,2𝑁−2), (28) where 𝐻⁽¹⁾_𝑁,𝑚≡𝑐_𝑚+1𝐻_𝑁,𝑚−𝑐_𝑚𝐻_𝑁,𝑚+1. (29) By this definition, 𝐻⁽¹⁾_𝑁,𝑚 are the Hankel matrices whose (1,1)-components are zero. We denote their components by 𝐻⁽¹⁾_𝑁,𝑚=0𝑐⁽¹⁾_𝑚+1⋯𝑐⁽¹⁾_𝑚+𝑁𝑐⁽¹⁾_𝑚+1𝑐⁽¹⁾_𝑚+2⋯𝑐⁽¹⁾_𝑚+𝑁+1𝑐⋮⋮⋮⋮_(1)𝑚+𝑁𝑐⁽¹⁾_𝑚+𝑁+1⋯𝑐⁽¹⁾_𝑚+2𝑁. (30) Next, from those 2𝑁−1 equations, we make 2𝑁−2 equations: 𝝈^𝑇𝐻⁽²⁾_𝑁,𝑚𝝈=0(𝑚=0,1,…,2𝑁−3), (31) where 𝐻⁽²⁾_𝑁,𝑚=𝑐⁽¹⁾_𝑚+2𝐻⁽¹⁾_𝑁,𝑚−𝑐⁽¹⁾_𝑚+1𝐻⁽¹⁾_𝑁,𝑚+1. (32) By definition, 𝐻⁽²⁾_𝑁,𝑚 are the Hankel matrices expressed by 𝐻⁽²⁾_𝑁,𝑚=00𝑐⁽²⁾_𝑚+2⋯𝑐⁽²⁾_𝑚+𝑁0𝑐⁽²⁾_𝑚+2⋯⋯𝑐⁽²⁾_𝑚+𝑁+1𝑐⁽²⁾_𝑚+2⋯⋯⋯𝑐⁽²⁾_𝑚+𝑁+1𝑐⋮⋮⋮⋮⋮_(2)𝑚+𝑁𝑐⁽²⁾_𝑚+𝑁+1⋯⋯𝑐⁽²⁾_𝑚+2𝑁. (33) Note that the components along the first and second lines in the antidiagonal direction become zero. Repeating this procedure to eliminate the components on the 𝑗th line in the anti-diagonal direction for 𝑗=3,…,2𝑁−1, we obtain 𝝈^𝑇𝐻^(𝑗)_𝑁,𝑚𝝈=0𝑚=0,1,…,2𝑁−1−𝑗, (34) where 𝐻^(𝑗)_𝑁,𝑚 are the Hankel matrices whose components along the first through 𝑗th lines in the anti-diagonal direction are zero. Here, (34) for 𝑗=2𝑁−1 is a linear equation for 𝜎_𝑁. Solving it and substituting the solution into (34) for 𝑗=2𝑁−2 gives a linear equation for 𝜎_𝑁−1. Repeating the backward substitution gives 𝜎_𝑁−2,…,𝜎₁, successively. The conditions for the unique solvability of these linear equations are 𝜁_𝑖≠𝜁_𝑗 for 𝑖≠𝑗 and 𝜈_𝑘≠0. See [17] for more details. Once 𝜎₁,…,𝜎_𝑁 are obtained, we have 𝜁_𝑘 as the roots of 𝜁^𝑁+𝜎₁𝜁^𝑁−1+⋯+𝜎_𝑁=0.

In a practical situation when data includes noise, the above procedure proposed in [17] to transform the second-degree equations (26) into linear equations (34) to obtain 𝜎_𝑁,…,𝜎₁ may be sensitive to the noise due to the cancellation. To reduce the cancellation error, we propose the following algorithm. First, we solve the second-degree equations (26) for 𝑚=0,…,𝑁−1 using a well-established method for solving simultaneous algebraic equations, for example, by means of the Gröbner bases. Then, out of the obtained 2^𝑁 solutions, (𝜎_1,𝑗,…,𝜎_𝑁,𝑗) for 𝑗=1,…,2^𝑁, we choose one which minimizes the sum of the absolute values of the left-hand sides of (26) for 𝑚=𝑁,…,2𝑁−1 given by ^2𝑁−1_𝑚=𝑁𝝈^𝑇𝐻_𝑁,𝑚𝝈. (35) In other words, to obtain the theoretically unique solution of (26) for 𝑚=0,…,2𝑁−1, first we find 2^𝑁 candidates from (26) for 𝑚=0,…,𝑁−1 and then choose the one which best explains the remaining equations (26) for 𝑚=𝑁,…,2𝑁−1.

3.2. Reconstruction of 𝜇_𝑘, 𝜈_𝑘, and 𝑁

To estimate 𝑁, following the method for the dipole source model [10], we assume that there are 𝑁^′(>𝑁) dipole-quadrupole sources and then estimate 𝜁_𝑘 for 𝑘=1,2,…,𝑁^′ using the algorithm in Section 3.1. Once they are obtained, 𝜇_𝑘 and 𝜈_𝑘 for 𝑘=1,2,…,𝑁^′ are linearly solved using (24) for 𝑛=0,1,…,2𝑁^′−1. Then, we compute |𝜇_𝑛+1/𝜇_𝑛| and |𝜈_𝑛+1/𝜈_𝑛| for 𝑛=1,2,…,𝑁^′−1, which are expected to be sufficiently small when 𝑛=𝑁. Practically, we estimate 𝑁 such that these ratios become smaller than some threshold set a priori. The thresholds should be determined by the ratios of the noise level contained in the data to the dipole and quadrupole strength which can be regarded as a physiologically meaningful source. As for the dipole source model in the 2D problem, the threshold is theoretically evaluated in the context of the Padé approximation [24]. A similar theory for the dipole-quadrupole source model, although greatly required, is left for further research; in this paper we show only numerical simulations in Section 4.

3.3. Algebraic Reconstruction of the 𝑧-Coordinates

To determine the 𝑧-coordinates of the dipole-quadrupole positions, after we determine the 𝑥𝑦-coordinates of the source, we use the method proposed in [7]. We put 𝑙=𝑛+1 and 𝑚=𝑛 in (16) and (17), where 𝑛∈ℕ^∗. From the identity 𝑟^𝑛+1𝑌_𝑛+1,𝑛=(2𝑛+1)!(𝑥+𝑖𝑦)^𝑛𝑧 [23], we have ^𝑁_𝑘=1∇𝑥+𝑖𝑦^𝑛𝑧_{𝐫=𝐫𝑘}⋅𝐫_𝑘×𝐩_𝑘+𝑄_𝑘∶∇∇𝑥+𝑖𝑦^𝑛𝑧×𝐫^𝑇_{𝐫=𝐫𝑘}=2𝑛+3𝜇_0Γ𝐧⋅𝐁(𝐫)𝑥+𝑖𝑦^𝑛𝑧𝑑𝑆. (36) Now, let 𝜙_𝑗(𝑡)≡𝑡−𝜁_𝑗^2𝑁_{𝑙=1,𝑙≠𝑗}𝑡−𝜁_𝑙³≡^3𝑁−1_𝑛=0𝛼_𝑛𝑡^𝑛. (37) Then, from ∑^3𝑁−1_𝑛=0𝛼_𝑛 × (18) and ∑^3𝑁−1_𝑛=0𝛼_𝑛 × (36), we obtain ^𝑁_𝑘=1∇𝜙_𝑗𝑥+𝑖𝑦_{𝐫=𝐫𝑘}⋅𝐫_𝑘×𝐩_𝑘+𝑄_𝑘∶∇∇𝜙_𝑗𝑥+𝑖𝑦×𝐫^𝑇_{𝐫=𝐫𝑘}=1𝜇₀^3𝑁−1_𝑛=0(2𝑛+1)𝛼_𝑛Γ𝐧⋅𝐁(𝐫)𝑥+𝑖𝑦^𝑛𝑑𝑆,_𝑁𝑘=1∇𝜙_𝑗𝑥+𝑖𝑦𝑧_{𝐫=𝐫𝑘}⋅𝐫_𝑘×𝐩_𝑘+𝑄_𝑘∶∇∇𝜙_𝑗𝑥+𝑖𝑦𝑧×𝐫^𝑇_{𝐫=𝐫𝑘}=1𝜇₀^3𝑁−1_𝑛=0(2𝑛+3)𝛼_𝑛Γ𝐧⋅𝐁(𝐫)𝑥+𝑖𝑦^𝑛𝑧𝑑𝑆, (38)(39) respectively. In the left-hand sides of (38) and (39), ∇𝜙_𝑗(𝑥+𝑖𝑦)_{𝐫=𝐫𝑘}=𝜙^′_𝑗𝜁_𝑘𝑖𝜙^′_𝑗𝜁_𝑘0=𝟎,∇𝜙_𝑗𝑥+𝑖𝑦𝑧_{𝐫=𝐫𝑘}=𝑧𝜙^′_𝑗𝜁_𝑘𝑖𝑧𝜙^′_𝑗𝜁_𝑘𝜙_𝑗𝜁_𝑘=𝟎 (40)(41) since from (37) it holds that 𝜙_𝑗(𝜁_𝑘)=𝜙^′_𝑗(𝜁_𝑘)=0 for all 𝑘=1,2,…,𝑁 irrespective of whether or not 𝑘 is equal to 𝑗. Also, we have ∇∇𝜙_𝑗𝑥+𝑖𝑦×𝐫_{𝐫=𝐫𝑘}=𝛿_𝑗𝑘Φ_𝑗,∇(∇𝜙_𝑗(𝑥+𝑖𝑦)𝑧×𝐫)_{𝐫=𝐫𝑘}=𝑧_𝑗𝛿_𝑗𝑘Φ_𝑗, (42)(43) where 𝛿_𝑗𝑘 is the Kronecker delta, and Φ_𝑗=2^𝑁_{𝑙=1,𝑙≠𝑗}𝜁_𝑗−𝜁_𝑙³𝑖𝑧_𝑗−𝑧_𝑗−𝑖𝜁_𝑗−𝑧_𝑗−𝑖𝑧_𝑗𝜁_𝑗000. (44) See Appendix D for the derivation of (42) and (43).

Substituting (40) through (43) into (38) and (39) gives 𝑄_𝑗∶Φ^𝑇_𝑗=1𝜇₀^3𝑁−1_𝑛=0(2𝑛+1)𝛼_𝑛Γ𝐧⋅𝐁(𝐫)𝑥+𝑖𝑦^𝑛𝑄𝑑𝑆,_𝑗∶𝑧_𝑗Φ^𝑇_𝑗=1𝜇₀^3𝑁−1_𝑛=0(2𝑛+3)𝛼_𝑛Γ𝐧⋅𝐁(𝐫)𝑥+𝑖𝑦^𝑛𝑧𝑑𝑆, (45) from which we obtain 𝑧_𝑗=∑^3𝑁−1_𝑛=0(2𝑛+3)𝛼_𝑛∫_Γ𝐧⋅𝐁(𝐫)𝑥+𝑖𝑦^𝑛𝑧𝑑𝑆∑^3𝑁−1_𝑛=0(2𝑛+1)𝛼_𝑛∫_Γ𝐧⋅𝐁(𝐫)𝑥+𝑖𝑦^𝑛𝑑𝑆. (46) This is the reconstruction formula of the 𝑧-coordinates: once 𝜁_𝑘 for 𝑘=1,2,…,𝑁 are obtained, then 𝑧_𝑗 are determined using (46) for each 𝑗=1,2,…,𝑁 with 𝛼_𝑛 defined by (37).

3.4. Algorithm

Our algorithm is summarized as follows.

Step 1 (estimate the number of the dipole-quadrupole sources). (1) Assume that there are 𝑁^′(>𝑁) dipole-quadrupole sources. Using the algorithm in Step 2 (1)–(5) below, where 𝑁 is replaced by 𝑁^′, obtain 𝜇_𝑘 and 𝜈_𝑘 for 𝑘=1,2,…,𝑁^′.
(2) Compute |𝜇_𝑛+1/𝜇_𝑛| and |𝜈_𝑛+1/𝜈_𝑛| for 𝑛=1,2,…,𝑁^′−1. When 𝑁|𝜇₊₁𝑁/𝜇|<𝜖_𝜇 and 𝑁|𝜈₊₁𝑁/𝜈|<𝜖_𝜈 for the thresholds 𝜖_𝜇 and 𝜖_𝜈, which are appropriately set a priori, we estimate 𝑁 to be 𝑁.

Step 2 (estimate the position parameters of 𝑁 dipole-quadrupole sources). (1) Given 𝐧⋅𝐁 on Γ, compute 𝑐_𝑛 in (23) for 𝑛=0,1,…,4𝑁−1 and construct the Hankel matrices (27) for 𝑚=0,1,…,2𝑁−1.
(2) Solve the simultaneous second-degree equation (26) for 𝑚=0,1,…,𝑁−1, using the Gröbner bases to get 2^𝑁 solutions.
(3) Out of the 2^𝑁 solutions, choose one solution 𝜎₁,…,𝜎_𝑁 that minimizes (35).
(4) Solve 𝑧^𝑁+𝜎₁𝑧^𝑁−1+⋯+𝜎_𝑁=0 to obtain 𝜁₁,…,𝜁_𝑁.
(5) Solve (24) for 𝜇_𝑘 and 𝜈_𝑘, where 𝑘=1,…,𝑁.
(6) Using the obtained 𝜁₁,…,𝜁_𝑁, compute 𝛼_𝑛 in (37) for 𝑛=0,1,…,3𝑁−1 and determine 𝑧₁,…,𝑧_𝑁 using (46).

4. Numerical Simulations

In this section, our algorithm is verified numerically. Γ was assumed to be a sphere with 𝑅=0.12m on which 361 measurement points were distributed uniformly using the spherical t-design [25]. To validate our algorithm for the dipole-quadrupole model, we assumed that the data was available on the whole sphere which enclosed the source. Identification using data on a part of Γ is an important aspect of future studies, for which the method proposed by Taulu et al. [22] would be useful.

4.1. Demonstration of Our Method

In this subsection, We examine the case where there are 𝑁=2 dipole-quadrupole sources whose parameters are given in Table 1. The theoretical data was computed using (11) to which Gaussian noise with 𝜎=10⁻² was added where the noise level 𝜎 is defined by the ratio of the standard deviation of the noise to the root mean squares of the data. 100 data sets with the different noise added were used for reconstruction.

Table 1 

Dipole-quadrupole parameters when 𝑁=2.

First, to determine the number of the sources 𝑁, we assume that there are 𝑁^′=3 dipole-quadrupole sources. Figure 2 shows the obtained |𝜇_𝑘| and |𝜈_𝑘| for 𝑘=1,2,3. One observes that |𝜇₃| and |𝜈₃| are much smaller than |𝜇_𝑘| and |𝜈_𝑘| for 𝑘=1,2, respectively. In fact, the geometric means of |𝜇₃/𝜇₂| and |𝜈₃/𝜈₂| for 100 data sets were 3.8×10⁻³ and 7.2×10⁻³, respectively. From this, we can judge that there are substantially 𝑁=2 poles.

(a)

(b)

Figure 2 

Estimation of 𝑁: the obtained |𝜇𝑘| (a) and |𝜈𝑘| (b) for 𝑘=1,2,3 when assuming 𝑁′=3. Since |𝜇3| and |𝜈3| are much smaller than |𝜇𝑘| and |𝜈𝑘| for 𝑘=1,2, respectively (|𝜇3/𝜇2|=3.8×10−3 and |𝜈3/𝜈2|=7.2×10−3), we can judge that there are substantially 𝑁=2 dipole-quadrupole sources.

Figure 3 shows the localization result projected on the 𝑥𝑦- and 𝑦𝑧-planes. We find that, under this noise level, two dipole-quadrupole sources were stably reconstructed. The average and standard deviations of the computational time to identify the parameters of the two dipole-quadrupole sources were 0.69 sec and 0.01 sec, respectively, using an Intel Core i7 CPU 870 (2.93 GHz) with 8 GB RAM, showing that our algebraic algorithm can determine the source parameters with low computational cost.

(a)

(b)

Figure 3 

Localization result when 𝑁=2 where the noise level 𝜎=10−2.

4.2. Effect of Noise on Localization Accuracy

To examine the stability of our algorithm, the noise level 𝜎 was varied in the wide range 10⁻⁸≤𝜎≤10⁻¹. The source and sensor configuration were the same as in Section 4.1. Figure 4 shows the standard deviation (s.d.) of the localization errors with respect to the noise level 𝜎 for the two poles. It is observed that s.d. increases as 𝜎 increases. For this source in the sphere with the radius 0.12 m, to realize s.d. less than 0.01 m, one finds that the noise level should be less than 𝜎=10⁻².

Figure 4 

Standard deviation of the localization errors with respect to the noise level 𝜎.

4.3. Effect of the Distance between Two Sources on Localization Accuracy

We varied the distance between the two dipole-quadrupole sources: the two dipole-quadrupole source were set at (0.02,±𝑑,0.07) where 2𝑑 was changed from 0.01 through 0.09. Figure 5 shows the relative localization error (error divided by the radius of Γ, 0.12). It is observed that the distance between the two dipole-quadrupole sources should be larger than about 0.06, that is the half of the radius of Γ, in order to obtain the relative error less than 10⁻².

Figure 5 

Relative localization error with respect to the distance between the two dipole-quadrupole sources.

4.4. Example When 𝑁=3

We examine the case when there are 𝑁=3 dipole-quadrupole sources whose parameters are given in Table 2. The noise level was set at 𝜎=10⁻². In this subsection, we assumed that 𝑁 is known a priori. Figure 6 shows the localization result. Comparing with the result when 𝑁=2 in Figure 3, the biases as well as the standard deviations of the localization errors become larger when 𝑁=3. It is numerically shown in [26] that, when using dipole-only source model, stability gets worse when 𝑁 increases and it is difficult to estimate the dipoles greater than about 𝑁=5. The example in this subsection shows that it is more difficult to estimate the large number (𝑁≥3) of dipole-quadrupole sources. This is considered due to fast spatial decay of the quadrupolar field. However, the dipole-quadrupole source even with 𝑁=1 or 2 has an advantage that it can well estimate a spatially distributed sources, especially when they are modeled with two, close, oppositely directed dipoles, which is a weak point of the conventional dipole model. This is shown in the next subsection.

Table 2 

dipole-quadrupole parameters when 𝑁=3.

(a)

(b)

Figure 6 

Localization result when 𝑁=3.

4.5. Identification of Dipoles Distributed on Curved Surfaces

In this subsection, we compare our algebraic method assuming the dipole-quadrupole model (DQM) with a conventional algebraic method assuming the dipole model (DM) for estimating the spatially distributed dipole sources. To model dipoles on cerebral convolutions, we assume that dipoles are placed on a mesh on a half-cylinder with a radius of 𝑟=5mm and a height of ℎ=5mm, as shown in Figure 7. There are six dipoles in the circumferential direction by five in the longitudinal direction, and hence a total of 30 dipoles on the half-cylinder. All the dipoles are aligned perpendicular to the surface of the cylinder to model the fact that the dipoles are perpendicular to the cerebral surface. We examined the following two cases.

Figure 7 

Source of distributed dipoles on a half-cylinder.

Case 1. A single half-cylinder source at 𝐫₁=(20,50,20) mm. The vectors which determine the posture of the cylinder, 𝐞₁₂ and 𝐞₁₁, are set to be 𝐞_𝑟 and 𝐞_𝜃, respectively, where 𝐞_𝑟 and 𝐞_𝜃 are the unit vectors in the 𝑟 and 𝜃 directions at 𝐫₁. (see Figure 8(a)). In this case, the total dipole moment 𝐩₁ is almost parallel to 𝐫₁; the angle between them is 2.4 degrees. Since a radial dipole is a silent source for the radial component of the magnetic field [1], this cylindrical source is regarded as being almost quadrupolar.

(a)

(b)

Figure 8 

Cases 1 and 2 of the cylindrically distributed dipoles.

Case 2. The same half-cylinder at 𝐫₁ as in Case 1, but with 𝐞₁₂=−𝐞_𝜃 and 𝐞₁₁=𝐞_𝑟 (Figure 8(b)). In this case, the angle between 𝐫₁ and 𝐩₁ is 87 degrees, so that the source has a detectable equivalent dipole moment as well as the equivalent quadrupole moment.

We computed the forward solution generated by the 30 elemental dipoles using (1). Note that (11) was not used to compute the theoretical data. 𝑁^′=3 was assumed to estimate 𝑁.

4.5.1. Case 1

Figure 9 shows the estimated |𝜇_𝑘| normalized by max{𝜇_𝑘} using DM. From the fact that |𝜇₃/𝜇₂|=0.005, 𝑁 is estimated to be two using DM. The reason why not a single equivalent dipole but two dipoles are estimated for a single cylindrical source is that the sum of the dipoles on the half-cylindrical surface is directed to the radial direction which is “silent” to the radial magnetic field outside the head. In contrast, Figures 10(a) and 10(b) show the estimated |𝜇_𝑘| normalized by max{𝜇_𝑘} and |𝜈_𝑘| normalized by max{𝜈_𝑘}, respectively, using DQM. From |𝜇₃/𝜇₂|=0.041 and |𝜈₃/𝜈₂|=0.004, 𝑁 is estimated to be one.

Figure 9 

Estimation of 𝑁 in Case 1 using the dipole model (DM) when assuming 𝑁′=3. From |𝜇3/𝜇2|=0.005, DM estimates 𝑁=2.

(a)

(b)

Figure 10 

Estimation of 𝑁 in Case 1 using the dipole-quadrupole model (DQM) when assuming 𝑁′=3. From |𝜇3/𝜇2|=0.041 and |𝜈3/𝜈2|=0.004, DQM estimates 𝑁=1.

Figure 11 shows the localization result using DM and DQM. It is observed that the two positions estimated with DM were far apart from the cylindrical surface, whereas DQM well estimated the center of the cylinder. This is a great advantage of DQM; although it is difficult for DM to estimate stably two, close, oppositely directed dipoles which are tangential to the spherical surface, it is easy for DQM.

(a)

(b)

Figure 11 

Localization result of the cylindrically distributed dipole sources using dipole model (DM) and dipole-quadrupole model (DQM).

4.5.2. Case 2

Figures 12, 13(a), and 13(b) show the estimated |𝜇_𝑘| using DM and |𝜇_𝑘| and |𝜈_𝑘| using DQM, respectively. In this Case, from |𝜇₂/𝜇₁|=0.006 in Figure 12 using DM and |𝜇₂/𝜇₁|=0.006 and |𝜈₂/𝜈₁|=0.073 in Figure 13 using DQM, both model estimate 𝑁=1. This is because the total of the elemental dipoles on the cylinder is almost perpendicular to the radial direction in Case 2 so that this cylindrical source can be regarded as a single equivalent current dipole. Figure 14 shows the localization result. Both DM and DQM can well estimate the center position of the cylindrical source in this case.

Figure 12 

Estimation of 𝑁 in Case 1 using the dipole model (DM) when assuming 𝑁′=3. From |𝜇3/𝜇2|=0.006, DM estimates 𝑁=1.

(a)

(b)

Figure 13 

Estimation of 𝑁 in Case 1 using the dipole-quadrupole model (DQM) when assuming 𝑁′=3. From |𝜇3/𝜇2|=0.006 and |𝜈3/𝜈2|=0.072, DQM estimates 𝑁=1.

(a)

(b)

Figure 14 

Localization result of the cylindrically distributed dipole sources using dipole model (DM) and dipole-quadrupole model (DQM).

5. Conclusion

We considered an inverse source problem of the Poisson equation for the radial component of the magnetic flux density when the current source is equivalently represented by multiple dipoles and quadrupoles in layered, spherical domains. By expressing the magnetic field in terms of the spherical harmonic expansion, we showed that the sectorial harmonics gave the simultaneous algebraic equations relating the dipole and quadrupole parameters to the boundary data. The equations had the same form as derived in the simple- and double-pole reconstruction problem in the 2D case [17], so that the source parameters can be algebraically reconstructed. We proposed a method to obtain the 𝑥- and 𝑦-coordinates of the dipole-quadrupole source by solving the simultaneous second-degree equations, with the result of which the 𝑧-coordinates are also determined. Numerical simulations show that the proposed algorithm estimates the number and centroid positions of the spatially extended dipoles well, especially when they include elemental dipoles oriented in opposite directions and their equivalent dipole moment is almost parallel to the radial direction of the spherical domain. A theoretical analysis to determine the thresholds depending on the noise in the algorithm for the estimation of the number of the dipole-quadrupole sources as well as an extension of the method to the case when the data is available only for part of the spherical boundary is left for further research.

Appendices

A. Derivation of (6)

Substituting the Taylor expansion of ∇^′(1/|𝐫−𝐫^′|) into the right-hand side of (3), we have (𝐫)=𝜇𝐫⋅𝐁⁰4𝜋^𝑁_{𝑘=1𝐷𝑘}∇^′1𝐫−𝐫^′_{𝐫′=𝐫𝑘}×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′+𝜇₀4𝜋^𝑁_{𝑘=1𝐷𝑘}∇^′∇^′1𝐫−𝐫^′_{𝐫′=𝐫𝑘}𝐫^′−𝐫_𝑘×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′+⋯. (A.1) The first term of (A.1) is rewritten as _𝐷𝑘∇^′1𝐫−𝐫^′_{𝐫′=𝐫𝑘}×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′=_𝐷𝑘𝐫−𝐫_𝑘𝐫−𝐫_𝑘³×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′=_𝐷𝑘𝐫−𝐫_𝑘𝐫−𝐫_𝑘³×𝐫_𝑘⋅𝐣_𝑝𝐫^′𝑑𝑣^′=𝐫_𝑘×_𝐷𝑘𝐣_𝑝𝐫^′𝑑𝑣^′⋅𝐫−𝐫_𝑘𝐫−𝐫_𝑘³. (A.2) The second term of (A.1) is rewritten as _𝐷𝑘∇^′∇^′1𝐫−𝐫^′_{𝐫′=𝐫𝑘}𝐫^′−𝐫_𝑘×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′=_𝐷𝑘∇^′𝐫−𝐫^′𝐫−𝐫^′3_{𝐫′=𝐫𝑘}𝐫^′−𝐫_𝑘×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′=_𝐷𝑘∇^′1𝐫−𝐫^′3𝐫−𝐫^′_{𝐫′=𝐫𝑘}𝐫^′−𝐫_𝑘×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′+_𝐷𝑘∇^′𝐫−𝐫^′𝐫−𝐫^′3_{𝐫′=𝐫𝑘}𝐫^′−𝐫_𝑘×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′=_𝐷𝑘3𝐫−𝐫^′𝐫−𝐫^′𝐫−𝐫^′5_{𝐫′=𝐫𝑘}𝐫^′−𝐫_𝑘×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′+_𝐷𝑘−𝐼𝐫−𝐫^′3_{𝐫′=𝐫𝑘}𝐫^′−𝐫_𝑘×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′=𝐼∶3×3identitymatrix_𝐷𝑘3𝐫−𝐫_𝑘𝐫−𝐫_𝑘⁵𝐫−𝐫_𝑘⋅𝐫^′−𝐫_𝑘×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′−_𝐷𝑘𝐫^′−𝐫_𝑘𝐫−𝐫_𝑘³×𝐫⋅𝐣_𝑝𝐫^′𝑑𝑣^′=_𝐷𝑘3𝐫×𝐫_𝑘𝐫−𝐫_𝑘𝐫−𝐫_𝑘⁵𝐫^′−𝐫_𝑘⋅𝐣_𝑝𝐫^′𝑑𝑣^′+_𝐷𝑘𝑋_𝑟𝐫^′−𝐫_𝑘𝐫−𝐫_𝑘³⋅𝐣_𝑝𝐫^′𝑑𝑣^′. (A.3) Here, for a 3×3 matrix 𝐴 whose row vectors are 𝐚^𝑇₁, 𝐚^𝑇₂, and 𝐚^𝑇₃, it holds that ()𝐚𝐴𝐱⋅𝐲=^𝑇1𝐱𝐚^𝑇₂𝐱𝐚^𝑇₃𝐱⋅𝐲=𝐚^𝑇₁𝐱𝑦₁+𝐚^𝑇₂𝐱𝑦₂+𝐚^𝑇₃𝐱𝑦₃=𝐴∶𝐲𝐱. (A.4) Using this property of the tensor contraction, (A.3) is further rewritten as 3𝐫×𝐫_𝑘𝐫−𝐫_𝑘𝐫−𝐫_𝑘⁵∶_𝐷𝑘𝐣_𝑝𝐫^′𝐫^′−𝐫_𝑘𝑑𝑣^′+𝑋_𝑟𝐫−𝐫_𝑘³∶_𝐷𝑘𝐣_𝑝𝐫^′𝐫^′−𝐫_𝑘𝑑𝑣^′=3𝐫×𝐫_𝑘𝐫−𝐫_𝑘+𝐫−𝐫_𝑘²𝑋_𝑟𝐫−𝐫_𝑘⁵∶_𝐷𝑘𝐣_𝑝𝐫^′𝐫^′−𝐫_𝑘𝑑𝑣^′. (A.5) Therefore, substituting (A.2) and (A.5) into (3) gives (6).

B.

This can be shown as follows. It is easy to obtain the dipole terms. For the quadrupole terms, we use the identity: _Ω1𝐚𝐫^′⋅𝑄_𝑘∇^′𝛿𝐫^′−𝐫_𝑘𝑑𝑣^′=𝑄∶∇^′𝐚𝐫^′𝑇_{𝐫′=𝐫𝑘} (B.1) for an arbitrary vector field 𝐚(𝐫^′)=(𝑎_𝑥(𝐫^′),𝑎_𝑦(𝐫^′),𝑎_𝑧(𝐫^′))^𝑇, where 𝑇 represents the transpose. When inserting the quadrupole terms in (12) into (1), we have from (B.1) _Ω1∇^′1𝐫−𝐫^′×𝐫^′⋅𝑄_𝑘∇^′𝛿𝐫^′−𝐫_𝑘𝑑𝑣^′=𝑄_𝑘∶∇^′∇^′1𝐫−𝐫^′×𝐫^′𝑇_{𝐫′=𝐫𝑘}. (B.2) Here, it holds that ∇^′∇^′1𝐫−𝐫^′×𝐫^′=∇^′𝐫×𝐫^′𝐫−𝐫^′3=∇^′1𝐫−𝐫^′3𝐫×𝐫^′+∇^′𝐫×𝐫^′𝐫−𝐫^′3=3𝐫−𝐫^′𝐫×𝐫^′𝐫−𝐫^′5+1𝐫−𝐫^′3,0𝑧−𝑦−𝑧0𝑥𝑦−𝑥0 (B.3) and hence ∇^′∇^′1𝐫−𝐫^′×𝐫^′𝑇_{𝐫′=𝐫𝑘}=3𝐫×𝐫_𝑘𝐫−𝐫_𝑘𝐫−𝐫_𝑘⁵+𝑋_𝑟𝐫−𝐫_𝑘³. (B.4) Thus we have (11).

C. Derivation of (21)

It holds that𝑄_𝑘∶∇∇𝑥+𝑖𝑦^𝑛×𝐫^𝑇_{𝐫=𝐫𝑘}=𝑞_{𝑥𝑥,𝑘}𝑞_{𝑥𝑦,𝑘}𝑞_{𝑥𝑧,𝑘}𝑞_{𝑥𝑦,𝑘}𝑞_{𝑦𝑦,𝑘}𝑞_{𝑦𝑧,𝑘}𝑞_{𝑥𝑧,𝑘}𝑞_{𝑦𝑧,𝑘}𝑞_{𝑧𝑧,𝑘}𝑖()𝑧∶𝑛𝑛−1𝑥+𝑖𝑦^𝑛−2−(𝑛−1)𝑧𝑥+𝑖𝑦^𝑛−2𝑖𝑥+𝑖𝑦^𝑛−1−(𝑛−1)𝑧𝑥+𝑖𝑦^𝑛−2−𝑖(𝑛−1)𝑧𝑥+𝑖𝑦^𝑛−2−𝑥+𝑖𝑦^𝑛−1−𝑖𝑛𝑥+𝑖𝑦^𝑛−1𝑛𝑥+𝑖𝑦^𝑛−10_{𝐫=𝐫𝑘}=𝑛𝑖(𝑛−1)𝑞_{𝑥𝑥,𝑘}𝑧_𝑘𝑥_𝑘+𝑖𝑦_𝑘^𝑛−2−(𝑛−1)𝑞_{𝑥𝑦,𝑘}𝑧_𝑘𝑥_𝑘+𝑖𝑦_𝑘^𝑛−2+𝑖𝑞_{𝑥𝑧,𝑘}𝑥_𝑘+𝑖𝑦_𝑘^𝑛−1−(𝑛−1)𝑞_{𝑥𝑦,𝑘}𝑧_𝑘𝑥_𝑘+𝑖𝑦_𝑘^𝑛−2−𝑖(𝑛−1)𝑞_{𝑦𝑦,𝑘}𝑧_𝑘𝑥_𝑘+𝑖𝑦_𝑘^𝑛−2−𝑞_{𝑦𝑧,𝑘}𝑥_𝑘+𝑖𝑦_𝑘^𝑛−1−𝑖𝑛𝑞_{𝑥𝑧,𝑘}𝑥_𝑘+𝑖𝑦_𝑘^𝑛−1+𝑛𝑞_{𝑦𝑧,𝑘}𝑥_𝑘+𝑖𝑦_𝑘^𝑛−1=𝑛(𝑛−1)𝑖𝑞_{𝑥𝑥,𝑘}−𝑞_{𝑦𝑦,𝑘}−2𝑞_{𝑥𝑦,𝑘}𝑧_𝑘+𝑞_{𝑦𝑧,𝑘}−𝑖𝑞_{𝑥𝑧,𝑘}𝑥_𝑘+𝑖𝑦_𝑘𝑥_𝑘+𝑖𝑦_𝑘^𝑛−2, (C.1)which gives (21).

D. Derivation of (42) and (43)

It is easy to see that ∇∇𝜙_𝑗𝑥+𝑖𝑦×𝐫_{𝐫=𝐫𝑘}=𝑖𝑧𝜙_𝑗^′′−𝑧𝜙_𝑗^′′−𝑖(𝑥+𝑖𝑦)𝜙_𝑗^′′−𝑖𝜙^′_𝑗−𝑧𝜙_𝑗^′′−𝑖𝑧𝜙_𝑗^′′(𝑥+𝑖𝑦)𝜙_𝑗^′′+𝜙^′_𝑗𝑖𝜙^′_𝑗−𝜙^′_𝑗0_{𝐫=𝐫𝑘}, (D.1) where 𝜙_𝑗, 𝜙^′_𝑗, and 𝜙_𝑗^′′ represent 𝜙_𝑗(𝑥+𝑖𝑦), 𝜙^′_𝑗(𝑥+𝑖𝑦), and 𝜙_𝑗^′′(𝑥+𝑖𝑦), respectively. From the definition of 𝜙_𝑗 in (37), one easily finds that 𝜙_𝑗(𝜁_𝑘)=𝜙^′_𝑗(𝜁_𝑘)=0 for all 𝑘=1,2,…,𝑁, and hence ∇∇𝜙_𝑗𝑥+𝑖𝑦×𝐫_{𝐫=𝐫𝑘}=𝑖𝑧_𝑘𝜙_𝑗^′′𝜁_𝑘−𝑧_𝑘𝜙_𝑗^′′𝜁_𝑘−𝑖𝜁_𝑘𝜙_𝑗^′′𝜁_𝑘−𝑧_𝑘𝜙_𝑗^′′𝜁_𝑘−𝑖𝑧_𝑘𝜙_𝑗^′′𝜁_𝑘𝜁_𝑘𝜙_𝑗^′′𝜁_𝑘.000 (D.2) Also, from the definition of 𝜙_𝑗, we have 𝜙_𝑗^′′𝜁_𝑘=2𝛿^𝑁_{𝑗𝑘𝑙=1,𝑙≠𝑗}𝜁_𝑗−𝜁_𝑙³, (D.3) where 𝛿_𝑗𝑘 is Kronecker’s delta, which gives (42). In the same way, ∇∇𝜙_𝑗𝑥+𝑖𝑦𝑧×𝐫_{𝐫=𝐫𝑘}=−𝑦𝜙^′_𝑗+𝑖𝑧²𝜙_𝑗^′′𝜙_𝑗+𝑥𝜙^′_𝑗−𝑧²𝜙_𝑗^′′−𝑖(𝑥+𝑖𝑦)𝑧𝜙_𝑗^′′−𝑖𝑧𝜙^′_𝑗−𝜙_𝑗−𝑖𝑦𝜙^′_𝑗−𝑧²𝜙_𝑗^′′𝑖𝑥𝜙^′_𝑗−𝑖𝑧²𝜙_𝑗^′′(𝑥+𝑖𝑦)𝑧𝜙_𝑗^′′+𝑧𝜙^′_𝑗2𝑖𝑧𝜙^′_𝑗−2𝑧𝜙^′_𝑗−𝑖(𝑥+𝑖𝑦)𝜙^′_{𝑗𝐫=𝐫𝑘}=𝑧_𝑘𝑖𝑧_𝑘𝜙_𝑗^′′𝜁_𝑘−𝑧_𝑘𝜙_𝑗^′′𝜁_𝑘−𝑖𝜁_𝑘𝜙_𝑗^′′𝜁_𝑘−𝑧_𝑘𝜙_𝑗^′′𝜁_𝑘−𝑖𝑧_𝑘𝜙_𝑗^′′𝜁_𝑘𝜁_𝑘𝜙_𝑗^′′𝜁_𝑘,000 (D.4)which gives (43).

Acknowledgments

This research has been partially supported by the Moritani Scholarship Foundation and the 2012 Inamori Research Grants Program.