Four-shell ellipsoidal model employing multipole expansion in ellipsoidal coordinates

Abstract

Although the head is more closely represented as an ellipsoid than a sphere, calculation in ellipsoidal coordinates is difficult. This paper presents a four shell ellipsoidal model, employing multipole expansion in ellipsoidal coordinates, for EEG, MEG, and evoked potential applications. Computational detail and insight into efficient calculation of the Lamé functions of the first and second kind are provided to demonstrate feasibilty. The Lamé function of the second kind, derived from the Lamé function of the first kind, can be computed at higher degrees by means of partial fraction expansion.

Appendices

Appendix A: Computation of external ellipsoidal harmonics using partial fraction expansion

The external ellipsoidal harmonics require the Lamé function of the second kind derived from Eq. (5). This equation can easily be determined for the first few degrees. However, as the degree, n, increases, the elliptic integral with the squared polynomial in the denominator becomes increasingly difficult to solve. To reduce the computational difficulty, the equation can be simplified by means of partial fraction expansion, each fraction containing an elliptic integral.

The roots of the Lamé functions of the first kind are obtained (e.g. using the Jenkins-Traub method [2, 32]), and the corresponding equation is expanded. The resulting partial fractions for each class of the Lamé functions are shown below. For odd n, K ^m _n (λ) is used to compute U ^m _n (λ)

$$ \begin{aligned} U_n^m(\lambda)&= (2n+1)K_n^m(\lambda) \int\limits_\lambda^\infty \frac{d\tau}{[\tau(\tau^n+a_{n-3}\tau^{n-2}+a_{n-5}\tau^{n-4} +\cdots+a_1)]^{2}\sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}}\\ &= (2n+1)K_n^m(\lambda)\left[ \int\limits_\lambda^\infty \frac{Ad\tau} {\tau^{2}\sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}} + \int\limits_\lambda^\infty \frac{Bd\tau}{(\tau^{2}-\gamma_1)^{2} \sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}}\right.\\ &\quad \left.+ \int\limits_\lambda^\infty \frac{Cd\tau}{(\tau^{2}-\gamma_2) \sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}} + \cdots \right] \end{aligned} $$

(22)

where A, B, C,... are the constants of the partial fractions, and γ’s are the roots from the K ^m _n (λ) polynomial. For n = 3, only the first three terms on the right side of Eq. (22) exist; for higher n, additional terms similar to those having B and C are present. For even n, Eq. (22) is the same except that A = 0.

Similarly, the expansion of V ^m _n using L ^m _n (λ) for even n is

$$ \begin{aligned} V_n^m(\lambda)&= (2n+1)L_n^m(\lambda)\int\limits_\lambda^\infty \frac{d\tau} {(\tau^{2}-b^{2}) [\tau(\tau^n+a_{n-2}\tau^{n-2}+a_{n-4}\tau^{n-4} +\cdots+a_0)]^{2}\sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}} \\& = (2n+1)L_n^m(\lambda)\left[\int\limits_\lambda^\infty \frac{Ad\tau} {\tau^{2}\sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}} + \int\limits_\lambda^\infty \frac{Bd\tau}{(\tau^{2}-\gamma_1)^{2} \sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}}\right.\\ &\quad \left.+ \int\limits_\lambda^\infty \frac{Cd\tau}{(\tau^{2}-\gamma_2) \sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}} + \cdots + \int\limits_\lambda^\infty {{Dd\tau}\over {(\tau^{2}-b^{2}) \sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}}}\right] \end{aligned} $$

(23)

The expansion of W ^m _n using M ^m _n (λ) is the same as Eq. (23) except that the terms involving (τ²−b ²) are replaced with those involving (τ²−c ²). For odd n, both V ^m _n (λ) and W ^m _n (λ) are the same for even n except that A = 0. The expansion of X ^m _n using N ^m _n (λ) for odd n is

$$ \begin{aligned} X_n^m(\lambda) = (2n+1)N_n^m\int\limits_\lambda^\infty \frac{d\tau} {(\tau^{2}-b^{2})(\tau^{2}-c^{2}) [\tau(\tau^n+a_{n-2}\tau^{n-2}+ \cdots+a_0)]^{2}\sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}}= (2n+1)N_n^m\left[\int\limits_\lambda^\infty \frac{Ad\tau} {\tau^{2}\sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}} + \int\limits_\lambda^\infty \frac{Bd\tau}{(\tau^{2}-\gamma_1)^{2} \sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}}\right. + \int\limits_\lambda^\infty \frac{Cd\tau}{(\tau^{2}-\gamma_2) \sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}} + \cdots + \int\limits_\lambda^\infty \frac{Dd\tau}{(\tau^{2}-b^{2}) \sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}}\left.\quad + \int\limits_\lambda^\infty \frac{Ed\tau}{(\tau^{2}-c^{2}) \sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}}\right] \end{aligned} $$

(24)

For even n, Eq. (24) can be used with A = 0.

Notice that the terms of all the classes (U ^m _n (x), V ^m _n (x), W ^m _n (x), and X ^m _n (x)) employ only five basic elliptic integrals for any degree. To simplify computations, these integrals can be replaced with the Carlson elliptic integrals [8] as follows:

$$ \begin{aligned} \int\limits_\lambda^\infty\frac{d\tau} {\tau^{2}\sqrt{\tau^{2}-b^{2}}\sqrt{\tau^{2}-c^{2}}}&= \frac{1} {3\lambda^3} R_D\left(1- \frac{u^{2}}{\lambda^{2}}, 1-\frac{v^{2}} {\lambda^{2}}, 1\right) \\ \int\limits_\lambda^\infty\frac{d\tau} {(\tau^{2}-b^{2})\sqrt{\tau^{2}-b^{2}} \sqrt{\tau^{2}-c^{2}}} &= \frac{1} {\lambda(b^{2}-c^{2})}\left[\sqrt{\frac{1-c^{2}/\lambda^{2}} {1-b^{2}/\lambda^{2}}}-R_F\left(1-\frac{u^{2}}{\lambda^{2}}, 1-\frac{v^{2}}{\lambda^{2}},1\right)\right. \\ &\quad \left.+\frac{c^{2}} {3\lambda^{2}}R_D\left(1-\frac{u^{2}}{\lambda^{2}}, 1-\frac{v^{2}} {\lambda^{2}},1\right)\right] \\ \int\limits_\lambda^\infty\frac{d\tau} {(\tau^{2}-c^{2})\sqrt{\tau^{2}-b^{2}} \sqrt{\tau^{2}-c^{2}}}= &\frac{1} {\lambda(c^{2}-b^{2})}\left[\sqrt{\frac{1-b^{2}/\lambda^{2}} {1-c^{2}/\lambda^{2}}}-R_F\left(1-\frac{u^{2}}{\lambda^{2}}, 1-\frac{v^{2}}{\lambda^{2}},1\right)\right. \\ &\quad \left.+\frac{b^{2}} {3\lambda^{2}}R_D\left(1-\frac{u^{2}}{\lambda^{2}}, 1-\frac{v^{2}} {\lambda^{2}},1\right)\right] \\ \int\limits_\lambda^\infty\frac{d\tau} {(\tau^{2}-\gamma)^{2}\sqrt{\tau^{2}-b^{2}} \sqrt{\tau^{2}-c^{2}}}= &\frac{1}{\lambda^3}\left\{\frac{1} {2(\gamma-b^{2})(\gamma-c^{2})} \left[\frac{\sqrt{\lambda^{2}-b^{2}}\sqrt{\lambda^{2}-c^{2}}} {1-\gamma/\lambda^{2}} \right.\right.\\ & \left.\quad+\frac{b^{2}c^{2}}{3\gamma} R_D\left(1-\frac{u^{2}}{\lambda^{2}},1-\frac{v^{2}} {\lambda^{2}},1\right)- \lambda^{2} R_F\left(1-\frac{u^{2}} {\lambda^{2}},1-\frac{v^{2}}{\lambda^{2}},1\right) \right] \\ &\quad \left.-\frac{1}{3\gamma}\left[1+\frac{\gamma^{2}-b^{2}c^{2}} {2(\gamma-b^{2})(\gamma-c^{2})}\right] R_J\left(1-\frac{u^{2}} {\lambda^{2}},1-\frac{v^{2}}{\lambda^{2}}, 1, 1-\frac{\gamma} {\lambda^{2}}\right)\right\} \\ \int\limits_\lambda^\infty\frac{d\tau} {(\tau^{2}-\gamma)\sqrt{\tau^{2}-b^{2}} \sqrt{\tau^{2}-c^{2}}}= &\frac{1}{3\lambda^3}R_J\left(1-\frac{u^{2}}{\lambda^{2}}, 1, 1-\frac{v^{2}}{\lambda^{2}},1-\frac{\gamma}{\lambda^{2}}\right) \end{aligned} $$

where R _F , R _D , and R _J are the Carlson elliptic integrals of the first, second, and third kinds, respectively. u = b and v = c if b > c, or u = c and v = b if c > b. Note that these equations require that b ≠ c. The first three integrals are for degree of one and higher, the fourth one for degree of two and higher, and the last one for degree of three and higher. Also, because each Lamé function of the second kind is derived from the corresponding Lamé function of the first kind, there are 2n + 1 functions of the second kind at each degree, n. Combining all these functions indicates that there are 4n + 2 functions of both kinds. To ensure that the internal and external ellipsoidal harmonics align at the boundaries, it is necessary to normalize the Lamé functions using the normalization integral developed in Appendix B.

Appendix B: Computation of normalization integral

To ensure that Eq. (15) holds, the functions E ^m _n (λ), E ^m _n (μ), E ^m _n (ν) must be normalized with the normalization integral [4, 28, 29]:

$$ c_{n,m}^4\int\limits_u^v\int\limits_0^u \frac{\left[\mu^{2}-\nu^{2}\right] \left[{\tilde{E}}_n^m(\mu){\tilde{E}}_n^m(\nu)\right]^{2} d\mu d\nu} {\sqrt{(\mu^{2}-b^{2}) (c^{2}-\mu^{2}) (b^{2}-\nu^{2}) (c^{2}-\nu^{2})}}=1 $$

(25)

where

$$ \begin{aligned} E_n^m(\lambda)&=c_{n,m}{\tilde{E}}_n^m(\lambda) \\ E_n^m(\mu)&=c_{n,m}{\tilde{E}}_n^m(\mu)\\ E_n^m(\nu)&=c_{n,m}{\tilde{E}}_n^m(\nu). \end{aligned} $$

(26)

${\tilde{E}}_n^m(\psi)$ and E ^m _n (ψ) are unnormalized and normalized, respectively. u = b and v = c if c > b or u = c and v = b if b > c. c _n,m are the normalization constants.

The normalization integral has the products of the functions ${\tilde{E}}_n^m(\mu){\tilde{E}}_n^m(\nu),$ which are equal to ${\tilde{K}}_n^m(\mu){\tilde{K}}_n^m(\nu), {\tilde{L}}_n^m(\mu){\tilde{L}}_n^m(\nu), {\tilde{K}}_n^m(\mu){\tilde{K}}_n^m(\nu),$ or ${\tilde{K}}_n^m(\mu){\tilde{K}}_n^m(\nu).$ Starting with the product

$$ {\tilde{E}}_n^m(\mu){\tilde{E}}_n^m(\nu)= {\tilde{K}}_n^m(\mu){\tilde{K}}_n^m(\nu), $$

the integral indicates that the product is squared. Squaring it results in

$$ \left[{\tilde{K}}_n^m(\mu){\tilde{K}}_n^m(\nu)\right]^{2}= \left\{\begin{array}{ll} \sum\limits_{\begin{subarray}{l}{i=0,1,2,\ldots,n}\\{j=0,1,2,\ldots, n}\end{subarray}} a_{ij}\mu^{2i}\nu^{2j} &\quad\hbox{for\ even}\ n\\ \sum\limits_{\begin{subarray}{l}{i=0,1,2,\ldots,n}\\{j=0,1,2,\ldots, n}\end{subarray}} a_{ij}\mu^{2(i+1)}\nu^{2(j+1)} &\quad\hbox{for\ odd}\ n \end{array}\right. $$

where a _ij ’s are the products of the coefficients of ${\tilde{K}}_n^m(\mu)\;\hbox{and}\;{\tilde{K}}_n^m(\nu).$ 2i, 2j, 2(i + 1), and 2(j + 1) are always even whether n is even or odd. These expressions are simplified because the squared product is intricate. Thus, each term in the normalization integral is

$$ \begin{aligned} \,&c_{n,m}^4\sum\limits_{\begin{subarray}{l}{i=0,1,2,\ldots,n}\\{j=0,1,2, \ldots,n}\end{subarray}}a_{ij}\int\limits_u^v\int\limits_0^u \frac{\mu^{2}-\nu^{2}}{\sqrt{(\mu^{2}-b^{2}) (c^{2}-\mu^{2}) (b^{2}-\nu^{2}) (c^{2}-\nu^{2})}} \mu^{2i}\nu^{2j} d\mu d\nu=1 \\ &\quad \quad \hbox{for\;even}\;n \\ &\, c_{n,m}^4 \sum\limits_{\begin{subarray}{l}{i=0,1,2,\ldots,n}\\{j=0,1,2,\ldots, n}\end{subarray}} a_{ij}\int\limits_u^v\int\limits_0^u \frac{\mu^{2}-\nu^{2}}{\sqrt{(\mu^{2}-b^{2}) (c^{2}-\mu^{2}) (b^{2}-\nu^{2}) (c^{2}-\nu^{2})}}\mu^{2(i+1)}\nu^{2(j+1)} d\mu d\nu=1 \\ &\quad \quad \hbox{for\;odd}\;n \\ \end{aligned} $$

where u = c and v = b if b > c, or u = b and v = c if c > b, and c _n,m is the normalization constant. Hence, the Lamé functions necessitate two normalization integrals: one is for odd n, and the other for even n.

Arranging the normalization integral,

$$ \begin{aligned} &c_{n,m}^4\sum\limits_{\begin{subarray}{l}{i=0,1,2,\ldots,n}\\{j=0,1,2, \ldots,n}\end{subarray}} a_{ij}[I_{i+1}G_j - I_iG_{j+1}]=1 \quad\hbox{for\;even}\;n\\ &c_{n,m}^4\sum\limits_{\begin{subarray}{l}{i=0,1,2,\ldots,n}\\{j=0,1,2, \ldots,n}\end{subarray}} a_{ij}[I_{i+2}G_{j+1} - I_{i+1}G_{j+2}]=1 \quad\hbox{for\;odd}\;n \end{aligned} $$

(27)

where

$$ \begin{aligned} I_i&=\int\limits_u^v \frac{\mu^{2i} d\mu} {\sqrt{|\mu^{2}-b^{2}|}\sqrt{|\mu^{2}-c^{2}|}} \quad\hbox{and}\\ G_j&=\int\limits_0^u\frac{\nu^{2j} d\nu} {\sqrt{b^{2}-\nu^{2}}\sqrt{{c^{2}-\nu^{2}}}}. \end{aligned} $$

Solving I and G with the elliptic integrals [7] and reducing them results in

$$ \left.\begin{array}{l}P_0=K(k^\prime)\\ P_1=E(k^\prime)\\ P_{i+1}=\frac{(1-2i)k^{2}P_{i-1}+2i(1+k^{2})P_i}{2i-1}\\ I_i=v^{2i-1}P_i \end{array} \right\} \quad \hbox{for the}\;I\;\hbox{integrals and} $$

(28)

$$ \left. \begin{array}{l}P_0=K(k)\\ P_1=K(k) - E(k)\\ P_{j+1}=\frac{(1-2j)k^{2}P_{j-1}+2j(1+k^{2})P_j}{2j-1}\\ G_j=v^{2j-1}P_j\end{array} \right\} \quad \hbox{for the}\;G\;\hbox{integrals} $$

(29)

where k = c/b if b > c or k = b/c if c > b, and $k^\prime=\sqrt{1-k^{2}}.$ K(k) and E(k) are the complete elliptic integrals of the first and second kinds, respectively. They can be replaced with the Carlson elliptic integrals [8]. Because P _i+1 and P _j+1 are recursive, I _i and G _j can easily be computed by programs at any degree n. Consequently, determining I _i and G _j yields c _n,m that normalize ${\tilde{K}}_n^m(\psi).$

Similarly,

$$ \left[{\tilde{L}}_n^m(\mu){\tilde{L}}_n^m(\nu)\right]^{2}= \sum\limits_{\begin{subarray}{l}{i=0,1,2,\ldots,n}\\{j=0,1,2,\ldots, n}\end{subarray}} a_{ij}|\mu^{2}-b^{2}||\nu^{2}-b^{2}|\mu^{2i}\nu^{2j} $$

where a _ij ’s are the products of the coefficients of ${\tilde{L}}_n^m(\mu)\;\hbox{and}\;{\tilde{L}}_n^m(\nu).$ For odd n, separating, arranging, and simplifying the equation leads to

$$ c_{n,m}^4\sum\limits_{\begin{subarray}{l}{i=0,1,2,\ldots,n}\\{j=0,1,2, \ldots,n}\end{subarray}} a_{ij}\left|(I_{i+1}G_{j+2}-I_{i+2}G_{j+1})+b^{2}(I_{i+2}G_j- I_iG_{j+2})+b^4(I_iG_{j+1}-I_{i+1}G_j)\right|=1. $$

(30)

For even n, the subscripts of I and G are incremented by one. Equation (30) is applicable to the functions ${\tilde{M}}_n^m(\mu){\tilde{M}}_n^m(\nu)$ except that b’s are replaced with c’s. Additionally, for ${\tilde{N}}_n^m(\mu){\tilde{N}}_n^m(\nu)$ of even n,

$$ \begin{aligned} \,&c_{n,m}^4\sum\limits_{\begin{subarray}{l}{i=0,1,2,\ldots,n} \\{j=0,1,2, \ldots,n}\end{subarray}} a_{ij}\left| I_{i+3}G_{j+2}-I_{i+2}G_{j+3}+ (b^{2}+c^{2})\left[I_{i+1}G_{j+3}-I_{i+3}G_{j+1}\right]\right.\\&\quad + b^{2}c^{2}[I_{i+2}G_{j+1}-I_{i+1}G_{j+2}+ I_{i+3}G_j-I_iG_{j+3}]+ b^{2}c^{2}(b^{2}+c^{2})[I_iG_{j+2}-I_{i+2}G_j]\\&\left.\quad + (b^4+c^4)[I_{i+2}G_{j+1}-I_{i+1}G_{j+2}]+ b^4c^4[I_{i+1}G_j-I_iG_{j+1}]\right|=1. \end{aligned} $$

Again, for odd n, the subscripts are incremented by one. Thus, solving I and G by means of P _i+1 and P _j+1 results in c _n,m that normalize the Lamé functions of the corresponding classes with Eq. (26).

In brief, the normalization integral normalizes the Lamé functions by setting their coefficients to obtain correct potentials. Arranging and simplifying the integral facilitates the computation of it by means of Eqs. (28) and (29) for computer programs. As a result, solving the normalization constant c _n,m leads to normalizing the functions with Eq. (26).