© The Authors 2026

1. Introduction

The magnetic field distribution in the deep interior of the Sun and the exact location of the solar dynamo remain active research topics (Charbonneau 2020; Fan 2021). These related features are driven by turbulent convective flows and large-scale dynamo mechanisms (Leighton et al. 1962; Weiss et al. 1996; Ossendrijver 2003; Miesch 2005; Fan 2009). The first significant changes in solar p-mode frequencies with the solar cycle were reported by Woodard & Noyes (1988).

Broomhall (2017) isolated the effect of the magnetic field on acoustic modes and showed that the frequencies of the solar p-modes vary in phase with the solar cycle by comparing observations between a maximum and a minimum of a solar cycle. Relatively few studies have attempted to use helioseismology to determine the internal magnetic field responsible for this frequency shift, although local magnetic features such as sunspots and near-surface flows are known to induce significant perturbations (Cally 2000; Schunker & Cally 2006; Gizon et al. 2009; Khomenko & Collados 2015; Rabello-Soares et al. 2018; Braun 2019).

General formalisms for computing multiplet frequency splittings were developed primarily to account for rotation and large-scale flow perturbations (Ritzwoller & Lavely 1991; Lavely & Ritzwoller 1992; Sekii et al. 1990; Sekii 1991; Schou et al. 1994; Hanasoge et al. 2017). Complementary local helioseismic techniques, such as ring-diagram analysis (Hill 1988; Basu et al. 1999) and time-distance helioseismology (Duvall et al. 1993; Giles 2000; Zhao & Kosovichev 2004; Hanasoge et al. 2012), have been used to infer subsurface flows and asphericities, providing independent constraints on the solar interior.

Early theoretical work on the impact of global magnetic fields on solar oscillations addressed the forward problem of frequency splitting due to rotation and axisymmetric internal magnetic fields (Gough & Thompson 1990; Dziembowski & Goode 1991). This approach was refined by Antia et al. (2000) to study magnetic fields in the convection zone. They set an upper limit of 300 kG for a toroidal magnetic field at the base of the convection zone. Baldner et al. (2009) extended this work to include a poloidal magnetic field component and compared the calculated splitting coefficients with observations over a solar cycle. They found that a combination of a poloidal and two shallow near-surface toroidal fields best reproduced the data.

Dziembowski & Goode (2004) further connected solar-cycle frequency shifts to near-surface turbulent velocity and small-scale magnetic fields, while also placing helioseismic limits on a buried toroidal field at the base of the convection zone. Antia et al. (2013) then used helioseismic measurements of angular-velocity variations to constrain the strength of internal magnetic fields. Hanasoge (2017) derived Lorentz-stress sensitivity kernels for general magnetic fields, and Das et al. (2020) further developed these results. Their framework, formulated within quasi-degenerate perturbation theory but reducible to the degenerate case for isolated multiplets, demonstrates the forward problem for both deep-toroidal and surface-dipolar magnetic configurations. These studies show that frequency shifts are particularly sensitive to near-surface perturbations and provide a formal basis for helioseismic constraints on internal magnetic fields.

Building on these helioseismic approaches, magnetoasteroseismology has explored how internal magnetic fields affect stellar oscillations beyond the Sun, particularly in red giants. In such stars a ‘fossil’ field can persist as a remnant of earlier dynamo processes (Braithwaite & Spruit 2004; Duez & Mathis 2010; Stello et al. 2016). Mathis & de Brye (2011) developed a formalism for magneto-gravito-inertial waves in magnetic, rotating stellar radiation zones, while Prat et al. (2019) applied similar ideas to compute gravity-mode frequency shifts in rapidly rotating stars. Several studies have applied these ideas to red giants, demonstrating that strong internal fields can alter mixed-mode frequencies (Bugnet et al. 2021), reduce asymptotic period spacings (Loi 2020), and depress dipole stellar oscillation modes (Fuller et al. 2015). More recent work has quantified core magnetic fields in red giants and explored the detectability of complex internal field topologies using mixed-mode asymmetries (Li et al. 2023; Das et al. 2024; Bhattacharya et al. 2024). Collectively, these works, which mostly focus on gravity modes and low-degree mixed modes (l = 1, 2), show that internal magnetic fields leave measurable signatures in stellar oscillations. This motivates extending helioseismic formalisms to asteroseismic contexts and complements our focus on p-mode perturbations.

In this context, Kiefer et al. (2017) and Kiefer & Roth (2018) provide a complementary theoretical approach and numerical results for the direct effect of a toroidal magnetic field on stellar oscillation frequencies, using quasi-degenerate perturbation theory (Lavely & Ritzwoller 1992) and properties of generalised spherical harmonics (GSHs, Dahlen & Tromp 1998). The direct effect is related to the inclusion of the Lorentz force in the wave equation, neglecting the indirect effect of the magnetic field on the stellar structure. We show that refining the assumptions used for calculating integrals over products of scalar GSHs reveals previously overlooked non-zero contributions, which can significantly affect the estimated frequency shift induced by the toroidal magnetic field.

While the approach of Das et al. (2020) applies to arbitrary fields, our work focuses on the fully analytical solution for a deep-seated toroidal field, providing explicit expressions and angular kernels that are particularly relevant to solar p-modes. In this paper, we refine the solution for the direct effect of a toroidal magnetic field on solar-like oscillations by presenting a thorough analytical method to compute integrals of scalar GSH products. We provide a corrected general matrix element and the corresponding angular kernels. Since a strong, large-scale toroidal magnetic field generated by a dynamo at the base of the convection zone is believed to produce active regions on the solar surface (Fan 2009; Charbonneau 2020), we recompute the shifts and splittings of multiplet frequencies induced by a strong, purely toroidal magnetic field at the base of the solar convection zone.

The structure of this paper is as follows. Section 2 describes the equilibrium stellar model perturbed by a magnetic field. Section 3 outlines the methods of quasi-degenerate perturbation theory. Section 4 presents the general matrix element for the direct effect of a toroidal magnetic field, the complete analytical solution for integrals over products of scalar GSHs, and the selection rules for mode coupling. Section 5 presents the shifts and splitting of unperturbed multiplet frequencies caused by a strong toroidal magnetic field at the base of the solar convection zone, which are discussed in Section 6.

2. Equilibrium model

In this section, we first introduce the eigenvalue problem for the oscillations of a non-magnetic, non-rotating star (Section 2.1). Following the approach of Kiefer et al. (2017), which is based on Lavely & Ritzwoller (1992), we then consider the magnetic field as a small perturbation to the oscillation equation (Section 2.2).

2.1. Non-magnetic, non-rotating equilibrium

The equation of motion for a gas cell in a non-rotating, non-magnetic star is given by Euler’s equation for an inviscid flow:

$$\begin{array}{c}\hfill \frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}·\mathrm{\nabla }\mathbf{v}=-\frac{1}{\rho }\mathrm{\nabla }p+\mathit{g},\end{array}$$(1)

where v is the local momentary velocity, p is the pressure, ρ is the density, and g is the gravitational acceleration. These reference stars are static and spherically symmetric. We present a brief derivation of the eigenvalue equation of adiabatic oscillations, whose solutions are the corresponding eigenfrequencies and eigenfunctions.

We introduced oscillations as small perturbations around the equilibrium state of the star, with equilibrium quantities labelled using the subscript 0. We assumed that all perturbations were sufficiently small to justify a linear approximation. We expressed the perturbations in the Eulerian form, denoted by a prime and defined for a quantity φ as

$$\begin{array}{c}\hfill \phi (\mathbf{r},t)={\phi }_0(\mathbf{r})+{\phi }^{\prime }(\mathbf{r},t).\end{array}$$(2)

The displacement of a fluid element is given by

$$\begin{array}{c}\hfill \mathit{\xi }=\mathbf{r}-{\mathbf{r}}_0,\end{array}$$(3)

and its velocity can be written as

$$\begin{array}{c}\hfill \mathbf{v}(\mathbf{r},t)=\frac{\partial \mathit{\xi }}{\partial t}.\end{array}$$(4)

By applying the linearisation given in Equation (2) to the equilibrium model’s quantities in the equation of motion (1) and using the velocity expression from Equation (4), we obtained

$$\begin{array}{c}\hfill {\rho }_0\frac{{\partial }^2\mathit{\xi }}{\partial t^2}={\rho }_0\frac{\partial \mathbf{v}}{\partial t}=-\mathrm{\nabla }{p}^{\prime }+{\rho }_0{\mathit{g}}^{\prime }+{\rho }^{\prime }{\mathit{g}}_0.\end{array}$$(5)

The gravitational acceleration, g, can be replaced with the gravitational potential, Φ, defined as g = −∇Φ, where Φ satisfies Poisson’s equation ∇²Φ = 4πGρ. Considering the spherical symmetry and time-independence of the equilibrium model, we adopted the following separation ansatz for the perturbation of a quantity φ:

$$\begin{array}{c}\hfill {\phi }^{\prime }(r,\theta ,\varphi ,t)={\stackrel{\sim }{\phi }}^{\prime }(r)Y_l^m(\theta ,\varphi ){\text{e}}^{-i\omega t},\end{array}$$(6)

where θ represents the co-latitude, ϕ denotes the longitude, Y_l^m(θ, ϕ) is a spherical harmonic characterised by its harmonic degree l and azimuthal order m, and ω is the angular frequency. We focused on a single eigenmode, characterised by its harmonic degree l, azimuthal order m, and eigenfrequency ω_nlm. Throughout the paper, the corresponding indices on perturbation quantities φ′ and ω are omitted for notational simplicity, unless needed for clarity. This separation removes the time dependence, allowing the linearised equation of motion to be expressed as an eigenvalue problem with eigenvalue ω² and a linear operator ℱ₀:

$$\begin{array}{c}\hfill {\omega }^2\mathit{\xi }=\frac{1}{{\rho }_0}\mathrm{\nabla }{p}^{\prime }-{\mathit{g}}^{\prime }-\frac{{\rho }^{\prime }}{{\rho }_0^2}\mathrm{\nabla }{p}_0={\mathcal{F}}_0(\mathit{\xi }),\end{array}$$(7)

where, in the non-magnetic case, we applied the equation of hydrostatic equilibrium ∇p₀ = ρ₀g₀. This form is convenient, as it generalises naturally to the equilibrium operator ${\stackrel{\sim }{\mathcal{F}}}_0$ in the presence of a magnetic field. This eigenvalue equation, combined with the appropriate boundary conditions, yields solutions only for specific discrete eigenvalues ω². For each (l, m), this provides a set of eigenfrequencies ω_nlm, characterised by their radial order n, which indicates the number of radial nodes (excluding a possible node at the centre) in the corresponding radial eigenfunction ${\stackrel{\sim }{\xi }}_{r,nlm}(r)$, as defined in Equation (11). At the surface, we impose zero-boundary conditions, assuming that both the pressure and the density vanish at r = R. Although a non-adiabatic atmosphere is generally necessary for completeness, it is excluded from this study for simplicity. The boundary conditions also eliminate singularities at the centre and surface, ensuring regularity. A detailed description of these boundary conditions can be found in Unno et al. (1989).

For a spherically symmetric star, the orientation of the coordinate system, which is tied to the azimuthal order m, has no physical meaning. Therefore, the eigenvalue ω_nlm² is (2l + 1)-fold degenerate with respect to m. The degenerate eigenfrequencies ω_nl obtained from the boundary value problem are referred to as multiplet frequencies. To solve the boundary eigenvalue problem and determine the eigenfrequencies and eigenfunctions, we first constructed a standard stellar model using the stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023). Using a MESA stellar evolution model as input, we used the pulsation code GYRE (Townsend & Teitler 2013; Townsend et al. 2018) to compute the eigenfrequencies and eigenfunctions for the adiabatic oscillation modes of the model. Instead of adopting Model S (Christensen-Dalsgaard et al. 1996), we constructed a standard solar model using MESA to illustrate that our method is applicable to general stellar models beyond the solar case. The root-mean-square of the relative difference in the sound speed between our MESA solar model and Model S is only 0.2%. The remaining discrepancies are primarily due to the omission of the atmospheric layers and a simplified treatment of material mixing. For a solar model, we computed modes up to 5200 μHz with harmonic degrees reaching l = 250.

In the following, we omit the tilde to distinguish between φ′(r, θ, ϕ, t) and $\stackrel{\sim }{\phi }\prime (r)$, as well as the subscript 0 for equilibrium quantities, for readability, unless needed for clarity. The displacement vector ξ, of which our eigenfunctions are composed, can be explicitly derived from Equation (7) using the separation of variables (6) and the gradient in spherical coordinates (A.4):

$$\begin{array}{cc}\hfill \mathit{\xi }=& {\xi }_r{\widehat{\mathbf{e}}}_r+{\mathit{\xi }}_h\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill =& \{\frac{1}{{\omega }^2}(\frac{1}{\rho }\frac{\partial p^{\prime }}{\partial r}+\frac{\partial {\mathrm{\Phi }}^{\prime }}{\partial r}+\frac{{\rho }^{\prime }}{\rho }g)Y_l^m{\widehat{\mathbf{e}}}_r\hfill \\ \hfill & +\frac{1}{r{\omega }^2}(\frac{p^{\prime }}{\rho }+{\mathrm{\Phi }}^{\prime })[\frac{\partial Y_l^m}{\partial \theta }{\widehat{\mathbf{e}}}_{\theta }+\frac{1}{sin\theta }\frac{\partial Y_l^m}{\partial \varphi }{\widehat{\mathbf{e}}}_{\varphi }]\}{\text{e}}^{-i\omega t}\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill =& \{{\stackrel{\sim }{\xi }}_r{Y}_l^m{\widehat{\mathbf{e}}}_r+{\stackrel{\sim }{\xi }}_h[\frac{\partial Y_l^m}{\partial \theta }{\widehat{\mathbf{e}}}_{\theta }+\frac{1}{sin\theta }\frac{\partial Y_l^m}{\partial \varphi }{\widehat{\mathbf{e}}}_{\varphi }]\}{\text{e}}^{-i\omega t}.\hfill \end{array}$$(10)

Here, we decompose the displacement vector into its radial and horizontal components, where ${\stackrel{\sim }{\xi }}_r(r)$ and ${\stackrel{\sim }{\xi }}_h(r)$ are defined as

$$\begin{array}{cc}\hfill {\stackrel{\sim }{\xi }}_r(r)=& \frac{1}{{\omega }^2}(\frac{1}{\rho }\frac{\partial p^{\prime }}{\partial r}+\frac{\partial {\mathrm{\Phi }}^{\prime }}{\partial r}+\frac{{\rho }^{\prime }}{\rho }g),\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {\stackrel{\sim }{\xi }}_h(r)=& \frac{1}{r{\omega }^2}(\frac{p^{\prime }}{\rho }+{\mathrm{\Phi }}^{\prime }).\hfill \end{array}$$(12)

2.2. Equilibrium model with magnetic field

The equation of motion (1) can be extended to include a magnetic field by adding the Lorentz force per mass (cgs units):

$$\begin{array}{c}\hfill \frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}·\mathrm{\nabla }\mathbf{v}=-\frac{1}{\rho }\mathrm{\nabla }p+\mathit{g}+\frac{1}{4\pi \rho }(\mathrm{\nabla }\times \mathbf{B})\times \mathbf{B},\end{array}$$(13)

where B is the magnetic field. We continue to assume a time-independent equilibrium model and neglect large-scale velocity fields, so that

$$\begin{array}{c}\hfill \frac{\partial {\stackrel{\sim }{\mathbf{B}}}_0}{\partial t}=0,{\stackrel{\sim }{\mathbf{v}}}_0(\mathbf{r})=0.\end{array}$$

To avoid confusion with the non-magnetic, non-rotating equilibrium quantities φ₀, we denote the equilibrium quantities in the presence of a magnetic field by ${\stackrel{\sim }{\phi }}_0$. Under equilibrium conditions, Equation (13) simplifies to

$$\begin{array}{c}\hfill 0=-\frac{1}{{\stackrel{\sim }{\rho }}_0}\mathrm{\nabla }{\stackrel{\sim }{p}}_0+{\stackrel{\sim }{\mathit{g}}}_0+\frac{1}{4\pi {\stackrel{\sim }{\rho }}_0}(\mathrm{\nabla }\times {\stackrel{\sim }{\mathbf{B}}}_0)\times {\stackrel{\sim }{\mathbf{B}}}_0.\end{array}$$(14)

We introduced stellar oscillations once more by linearising all quantities φ around the equilibrium model in the presence of a magnetic field,

$$\begin{array}{c}\hfill \phi (\mathbf{r},t)={\stackrel{\sim }{\phi }}_0(\mathbf{r})+{\phi }^{\prime }(\mathbf{r},t).\end{array}$$(15)

By applying linearisation to Equation (13), retaining only linear-order perturbations, separating the time dependencies as in Equation (6), and incorporating Equation (14), we obtained

$$\begin{array}{c}\hfill {\omega }^2\mathit{\xi }={\stackrel{\sim }{\mathcal{F}}}_0(\mathit{\xi })+{\mathcal{F}}_B(\mathit{\xi }).\end{array}$$(16)

The operator ${\stackrel{\sim }{\mathcal{F}}}_0$ is defined analogously to ℱ₀ given in Equation (7), but with equilibrium quantities ${\stackrel{\sim }{\phi }}_0$. The operator ℱ_B represents the direct effect of the magnetic field and is given by

$$\begin{array}{cc}\hfill {\mathcal{F}}_B(\mathit{\xi }):=& -\frac{1}{4\pi {\stackrel{\sim }{\rho }}_0}[\{\frac{(\mathit{\xi }·\mathrm{\nabla }){\stackrel{\sim }{\rho }}_0}{{\stackrel{\sim }{\rho }}_0}+\mathrm{\nabla }·\mathit{\xi }\}(\mathrm{\nabla }\times {\stackrel{\sim }{\mathbf{B}}}_0)\times {\stackrel{\sim }{\mathbf{B}}}_0\hfill \\ \hfill & +(\mathrm{\nabla }\times {\stackrel{\sim }{\mathbf{B}}}_0)\times {\mathbf{B}}^{\prime }+(\mathrm{\nabla }\times {\mathbf{B}}^{\prime })\times {\stackrel{\sim }{\mathbf{B}}}_0],\hfill \end{array}$$(17)

where the Eulerian perturbation of the magnetic field B′ is expressed using ideal magnetohydrodynamics as

$$\begin{array}{c}\hfill {\mathbf{B}}^{\prime }=\mathrm{\nabla }\times (\mathit{\xi }\times {\stackrel{\sim }{\mathbf{B}}}_0).\end{array}$$(18)

This result is consistent with Equation (19.13) in Unno et al. (1989) and Equation (2.5) in Gough & Thompson (1990). However, our result differs from Equation (24) in Kiefer et al. (2017). Kiefer et al. (2017) appears to have cancelled the density ${\stackrel{\sim }{\rho }}_0$ in the first term within the curly brackets of Equation (17). This approach is incorrect because the gradient operator in the numerator still acts on the density. Eliminating the density neglects the effect of the gradient on ${\stackrel{\sim }{\rho }}_0$, which is crucial for the correct formulation of the operator. Consequently, the analytical expressions for the operator ℱ_B and the general matrix element defined in Section 4 differ from those of Kiefer et al. (2017).

Equations (14) and (16) show that the magnetic field affects both the equilibrium structure and the oscillations. We considered both effects as small perturbations. To account for the change in the equilibrium structure we expanded the equilibrium quantities in the presence of a magnetic field, ${\stackrel{\sim }{\phi }}_0$, to first order around the non-magnetic, non-rotating equilibrium model given by the operator ℱ₀, as follows:

$$\begin{array}{c}\hfill {\stackrel{\sim }{\phi }}_0(r,\theta ,\varphi )={\phi }_0(r)+{\phi }_1(r,\theta ,\varphi ).\end{array}$$(19)

This gives

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathcal{F}}}_0(\mathit{\xi })={\mathcal{F}}_0(\mathit{\xi })+{\mathcal{F}}_1(\mathit{\xi }),\end{array}$$(20)

where ℱ₁ contains all linear-order perturbations to the stellar structure due to the magnetic field and is referred to as the indirect effect. The direct contribution of the magnetic field perturbation is given by ℱ_B, which couples the modes and thus modifies the oscillation frequencies and eigenfunctions. Although ℱ_B includes equilibrium structure quantities in the presence of the magnetic field, we retained only the leading-order term, since ℱ_B is also considered a small perturbation to the oscillation equations. Finally, we obtained the total adiabatic oscillation eigenvalue equation

$$\begin{array}{c}\hfill {\omega }^2\mathit{\xi }={\mathcal{F}}_0(\mathit{\xi })+{\mathcal{F}}_1(\mathit{\xi })+{\mathcal{F}}_B(\mathit{\xi })=\mathcal{F}(\mathit{\xi }).\end{array}$$(21)

The direct and indirect effects can be separated due to the additivity of the contributions, as shown in Equation (21). Kiefer & Roth (2018) provide a calculation of the indirect effect due to a magnetic field. In this work, we focus on the direct effect of the magnetic field ℱ_B. We therefore neglected the indirect effect ℱ₁ in Equation (21) and assumed that the equilibrium state is described by the non-magnetic, non-rotating model given by

$$\begin{array}{c}\hfill {\omega }_{0,k}^2{\mathit{\xi }}_{0,k}={\mathcal{F}}_0({\mathit{\xi }}_{0,k}),\end{array}$$(22)

where k = (n, l, m).

3. Quasi-degenerate perturbation theory

The eigenfrequencies of the equilibrium model form a dense spectrum, resulting in small frequency differences between the modes relative to the oscillation frequency. Therefore, we can apply the perturbation theory of quasi-degenerate states (Lavely & Ritzwoller 1992; Roth 2002; Kiefer et al. 2017). The eigenspace K consists of eigenfunctions with nearly identical eigenfrequencies. The eigenspace K includes an eigenfunction ξ_0, k only if its degenerate eigenfrequency ω_0, k² = ω_0, nl² satisfies

$$\begin{array}{c}\hfill |{\omega }_{0,k}^2-{\omega }_{\mathrm{ref}}^2|<\mathrm{\Delta }{\omega }^2\text{for}k\in K,\end{array}$$(23)

where ω_ref² is a reference frequency set equal to the degenerate frequency ω_0, j² of the multiplet that we want to perturb. The frequency range of the coupling modes within the eigenspace K is determined by Δω² and provides a measure of the accuracy of the perturbed quantities.

We assumed that the change in the oscillation equation due to the direct effect of the magnetic field was a small perturbation. We therefore expanded the eigenfunctions and eigenfrequencies of a mode j to first-order as

$$\begin{array}{cc}& \mathcal{F}={\mathcal{F}}_0+{\mathcal{F}}_B,\hfill \end{array}$$(24)

$$\begin{array}{cc}& {\mathit{\xi }}_j={\stackrel{\sim }{\mathit{\xi }}}_{0,j}+{\stackrel{\sim }{\mathit{\xi }}}_{1,j},\hfill \end{array}$$(25)

$$\begin{array}{cc}& {\omega }_j^2={\omega }_{0,j}^2+{\omega }_{1,j}^2.\hfill \end{array}$$(26)

The subscripts 0 and 1 in the expansion terms refer to the zeroth- and first-order perturbations, respectively. Here, ω_1, j² represents the first-order perturbation of the squared eigenfrequency of mode j. We further expanded the zeroth- and first-order eigenfunction perturbations in terms of the unperturbed eigenfunctions according to the quasi-degenerate perturbation theory as follows:

$$\begin{array}{cc}& {\stackrel{\sim }{\mathit{\xi }}}_{0,j}={\displaystyle \sum _{k\in K}}{a}_{\mathit{jk}}{\mathit{\xi }}_{0,k},\hfill \end{array}$$(27)

$$\begin{array}{cc}& {\stackrel{\sim }{\mathit{\xi }}}_{1,j}={\displaystyle \sum _{k\in K^{\perp }}}{b}_{\mathit{jk}}{\mathit{\xi }}_{0,k},\hfill \end{array}$$(28)

where K^⊥ is the orthogonal complement of K in the set of all unperturbed eigenfunctions. The expansion coefficients a_jk and b_jk are complex constants. We inserted the expansions (24)–(28) into the not-exactly solvable oscillation eigenvalue problem (21), discarding second- or higher-order perturbations and neglecting the indirect effect. We then computed the inner product with ρ₀ξ_0, k′^* and the resulting expression is integrated over the volume of the star. This yields

$$\begin{array}{cc}\hfill 0=& {\displaystyle \sum _{k\in K}}{a}_{\mathit{jk}}({\displaystyle {\int }_V}{\mathit{\xi }}_{0,k^{\prime }}^{\ast }({\mathcal{F}}_0-{\omega }_{0,j}^2){\mathit{\xi }}_{0,k}{\rho }_0\text{d}V+{\displaystyle {\int }_V}{\mathit{\xi }}_{0,k^{\prime }}^{\ast }({\mathcal{F}}_B-{\omega }_{1,j}^2)\hfill \\ \hfill & ·{\mathit{\xi }}_{0,k}{\rho }_0\text{d}V)+{\displaystyle \sum _{k\in K^{\perp }}}{b}_{\mathit{jk}}{\displaystyle {\int }_V}{\mathit{\xi }}_{0,k^{\prime }}^{\ast }({\mathcal{F}}_0-{\omega }_{0,j}^2){\mathit{\xi }}_{0,k}{\rho }_0\text{d}V.\hfill \end{array}$$(29)

The unperturbed eigenfunctions satisfy the orthogonality relation (Unno et al. 1989; Lavely & Ritzwoller 1992)

$$\begin{array}{c}\hfill {\displaystyle {\int }_V}{\rho }_0{\mathit{\xi }}_{k^{\prime }}^{\ast }·{\mathit{\xi }}_k\text{d}V={\delta }_{k^{\prime }k}{N}_k,\end{array}$$(30)

where we used the following abbreviations:

$$\begin{array}{cc}& N_k={\displaystyle {\int }_0^R}{\rho }_0({\xi }_r^{{}^{\prime }\ast }{\xi }_r+l(l+1){\xi }_h^{{}^{\prime }\ast }{\xi }_h)r^2\text{d}r,\hfill \end{array}$$(31)

$$\begin{array}{cc}& {\delta }_{k^{\prime }k}={\delta }_{l^{\prime }l}{\delta }_{m^{\prime }m}{\delta }_{n^{\prime }n}.\hfill \end{array}$$(32)

Using this orthogonality, together with Equation (22), and incorporating certain factors into the expansion coefficients, we derive

$$\begin{array}{c}\hfill {\displaystyle \sum _{k\in K}}{a}_{\mathit{jk}}{\omega }_{1,j}^2{\delta }_{k^{\prime }k}={\displaystyle \sum _{k\in K}}{a}_{\mathit{jk}}{Z}_{k^{\prime }k}-{\displaystyle \sum _{k\in K^{\perp }}}{b}_{\mathit{jk}}({\omega }_{0,j}^2-{\omega }_{0,k}^2){\delta }_{k^{\prime }k},\end{array}$$(33)

where Z_k′k is a component of the so-called supermatrix Z, defined for (k, k′) ∈ K as

$$\begin{array}{c}\hfill {\displaystyle Z_{k^{\prime }k}=\frac{1}{N_k}({\displaystyle {\int }_V}{\mathit{\xi }}_{0,k^{\prime }}^{\ast }{\mathcal{F}}_B({\mathit{\xi }}_{0,k}){\rho }_0\text{d}V-({\omega }_{0,j}^2-{\omega }_{0,k}^2)N_k{\delta }_{k^{\prime }k}).}\end{array}$$(34)

The last term in Equation (33) is non-zero only when k = k′∈K^⊥, which decouples it from the remaining terms for k ∈ K. This yields a decoupled eigenvalue equation,

$$\begin{array}{c}\hfill {\displaystyle \sum _{k\in K}}{a}_{\mathit{jk}}{Z}_{k^{\prime }k}={\displaystyle \sum _{k\in K}}{a}_{\mathit{jk}}{\omega }_{1,j}^2{\delta }_{k^{\prime }k},\end{array}$$(35)

which can also be expressed as

$$\begin{array}{c}\hfill \mathit{Z}\mathit{A}=\mathit{A}\mathbf{\Lambda },\end{array}$$(36)

where Λ is a diagonal matrix and A is the eigenvector matrix. The entries of Λ are the squared eigenfrequency perturbations, Λ_kj = ω_1, j²δ_kj for j ∈ K. Each column of A corresponds to an eigenvector a_j that contains the expansion coefficients such that A_kj = a_jk. Numerical methods can solve the eigenvalue Equation (36) for a given supermatrix Z. We define the general matrix element based on Equation (34) as

$$\begin{array}{c}\hfill H_{n\prime n,l\prime l}^{m\prime m}\equiv H_{k\prime k}={\displaystyle {\int }_V}{\mathit{\xi }}_{0,k\prime }^{\ast }{\mathcal{F}}_B({\mathit{\xi }}_{0,k}){\rho }_0\text{d}V.\end{array}$$(37)

This element represents the (m′,m) component of the general matrix H_{n′n, l′l}, where |m′| ≤ l′ and |m|≤l. Therefore, the component Z_k′k of the supermatrix Z can also be written as

$$\begin{array}{c}\hfill Z_{k^{\prime }k}=\frac{1}{N_k}(H_{k\prime k}-({\omega }_{0,j}^2-{\omega }_{0,k}^2)N_k{\delta }_{k\prime k}).\end{array}$$(38)

Solving the eigenvalue Equation (36) numerically yields the squared frequency perturbations ω_1, j² and the eigenvectors a_j. We used these quantities to determine the zeroth-order perturbation of the eigenfunctions ${\stackrel{\sim }{\mathit{\xi }}}_{0,j}$ from Equation (27). We matched the frequency shift ω_1, j² with the corresponding quantum numbers (n, l, m) of the unperturbed mode by extracting the maximum component a_j, kmax of the eigenvector a_j. Since the perturbations are small, we assumed that the main contribution to the perturbed eigenfunction ${\stackrel{\sim }{\mathit{\xi }}}_{0,j}$ originates from the unperturbed eigenfunction with the same quantum numbers. Thus, the degenerate eigenfrequency ω_0, x² corresponding to the mode ξ_0, kmax is subject to the frequency shift ω_1, j². The perturbed squared eigenfrequencies are then given by

$$\begin{array}{c}\hfill {\omega }_j^2={\omega }_{0,x}^2+{\omega }_{1,j}^2.\end{array}$$(39)

The frequency shift δω_j to a mode j is defined as (Lavely & Ritzwoller 1992):

$$\begin{array}{c}\hfill \delta {\omega }_j={\omega }_j-{\omega }_{0,x}.\end{array}$$(40)

Substituting this into Equation (39) and assuming a small perturbation while neglecting higher-order terms yields the following approximation for the frequency shift:

$$\begin{array}{c}\hfill \delta {\omega }_j=\frac{{\omega }_{1,j}^2}{2{\omega }_{0,x}}.\end{array}$$(41)

4. General matrix element for a toroidal magnetic field

We derived the general matrix element corresponding to a perturbation induced by a toroidal magnetic field following the approach outlined by Kiefer et al. (2017). Our analytical results differ from those presented by Kiefer et al. (2017) because we derived a different perturbation operator (Equation (17)). We obtained the general matrix element for a magnetic field by combining Equations (17) and (37):

$$\begin{array}{cc}\hfill H_{k\prime k}=& -\frac{1}{4\pi }{\displaystyle {\int }_V}{\mathit{\xi }}_{0,k\prime }^{\ast }[\{\frac{({\mathit{\xi }}_{0,k}·\mathrm{\nabla }){\rho }_0}{{\rho }_0}+\mathrm{\nabla }·{\mathit{\xi }}_{0,k}\}(\mathrm{\nabla }\times {\mathbf{B}}_0)\times {\mathbf{B}}_0\hfill \\ \hfill & +(\mathrm{\nabla }\times {\mathbf{B}}_0)\times {\mathbf{B}}^{\prime }+(\mathrm{\nabla }\times {\mathbf{B}}^{\prime })\times {\mathbf{B}}_0]\text{d}V.\hfill \end{array}$$(42)

Here, we omit the tilde denoting the equilibrium model under the influence of the magnetic field to focus on the direct effect and avoid confusion with the later use of the tilde notation. We modelled the magnetic field as a superposition of purely toroidal fields without longitudinal dependence. In practice, we achieved this by setting m = 0 in the generalised spherical harmonics (GSHs) of the magnetic field. Appendix D introduces the scalar GSHs. Thus, the toroidal magnetic field is expressed as

$$\begin{array}{c}\hfill \mathbf{B}={\mathbf{B}}_{\mathrm{tor}}={\displaystyle \sum _s}{\mathbf{B}}_s,\end{array}$$(43)

where

$$\begin{array}{c}\hfill {\mathbf{B}}_s=-a_s(r)\frac{\partial }{\partial \theta }{Y}_s^{0,0}(\theta ,\varphi ){\widehat{\mathbf{e}}}_{\varphi }.\end{array}$$(44)

The radial profile of the toroidal field components with degree s is denoted by a_s(r), while the angular profile is given by the co-latitudinal derivative of the scalar GSH $Y_s^{0,0}(\theta ,\varphi )$. The sum over the harmonic degrees s includes only even integers to ensure the observed opposite polarities in the northern and southern hemispheres (Hathaway 2015). By inserting the toroidal magnetic field into Equation (42), applying the distributivity of vector and dot products, using the sum rule for integrals, and replacing the Eulerian perturbation of the magnetic field with Equation (18), we obtained the following:

$$\begin{array}{cc}\hfill H_{k\prime k}=& -\frac{1}{4\pi }{\displaystyle \sum _{s,s\prime }}{\displaystyle {\int }_V}{\mathit{\xi }}_{0,k\prime }^{\ast }[\{\frac{({\mathit{\xi }}_{0,k}·\mathrm{\nabla }){\rho }_0}{{\rho }_0}+\mathrm{\nabla }·{\mathit{\xi }}_{0,k}\}(\mathrm{\nabla }\times {\mathbf{B}}_s)\times {\mathbf{B}}_{s\prime }\hfill \\ \hfill & +(\mathrm{\nabla }\times {\mathbf{B}}_s)\times \{\mathrm{\nabla }\times ({\mathit{\xi }}_{0,k}\times {\mathbf{B}}_{s\prime })\}\hfill \\ \hfill & +(\mathrm{\nabla }\times \{\mathrm{\nabla }\times ({\mathit{\xi }}_{0,k}\times {\mathbf{B}}_s)\})\times {\mathbf{B}}_{s\prime }]\text{d}V,\hfill \end{array}$$(45)

where s and s′ span the same index range.

We now focus on a single summand with fixed values of s and s′ and omit the indices s and s′ using the following notations:

$$\begin{array}{cc}\hfill {\mathbf{B}}_s& =-B_{\varphi }{\widehat{\mathbf{e}}}_{\varphi }=-a(r)\frac{\partial }{\partial \theta }{Y}_s^{0,0}(\theta ,\varphi ){\widehat{\mathbf{e}}}_{\varphi },\hfill \end{array}$$(46)

$$\begin{array}{cc}\hfill {\mathbf{B}}_{s\prime }& =-{\stackrel{\sim }{B}}_{\varphi }{\widehat{\mathbf{e}}}_{\varphi }=-\stackrel{\sim }{a}(r)\frac{\partial }{\partial \theta }{Y}_{s\prime }^{0,0}(\theta ,\varphi ){\widehat{\mathbf{e}}}_{\varphi },\hfill \end{array}$$(47)

$$\begin{array}{cc}\hfill {\mathit{\xi }}_{0,k}& ={\xi }_r(r)Y_l^{0,m}{\widehat{\mathbf{e}}}_r+{\xi }_h(r)[\frac{\partial }{\partial \theta }{Y}_l^{0,m}{\widehat{\mathbf{e}}}_{\theta }+\frac{1}{sin\theta }\frac{\partial }{\partial \varphi }{Y}_l^{0,m}{\widehat{\mathbf{e}}}_{\varphi }],\hfill \end{array}$$(48)

$$\begin{array}{cc}\hfill {\mathit{\xi }}_{0,k\prime }^{\ast }& ={\stackrel{\sim }{\xi }}_r(r)Y_{l\prime }^{0,m\prime \ast }{\widehat{\mathbf{e}}}_r+{\stackrel{\sim }{\xi }}_h(r)[\frac{\partial }{\partial \theta }{Y}_{l\prime }^{0,m\prime \ast }{\widehat{\mathbf{e}}}_{\theta }+\frac{1}{sin\theta }\frac{\partial }{\partial \varphi }{Y}_{l\prime }^{0,m\prime \ast }{\widehat{\mathbf{e}}}_{\varphi }].\hfill \end{array}$$(49)

We chose ξ_r and ξ_h to be real-valued (Aerts et al. 2010). A straightforward but lengthy calculation provides an explicit expression for the operator of the direct effect of a toroidal magnetic field perturbation ℱ_B(ξ_0, k), for fixed values of s and s′. Appendix B presents the operator for the direct effect of the toroidal magnetic field and the necessary computations to obtain the corresponding general matrix element. The explicit expression of the general matrix element for a fixed s and s′ is given by

$$\begin{array}{cc}\hfill H_{k\prime k}=& {\displaystyle {\int }_V}{\mathit{\xi }}_{0,k\prime }^{\ast }{\mathcal{F}}_B({\mathit{\xi }}_{0,k}){\rho }_0\text{d}V\hfill \\ \hfill =& \frac{1}{4\pi }{\displaystyle {\int }_V}[(\frac{a\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_r}{r}\frac{\partial ln{\rho }_0}{\partial r}+\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_r\frac{\partial ln{\rho }_0}{\partial r}\frac{\partial a}{\partial r}+\frac{a\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_r}{r^2}\hfill \\ \hfill & -\frac{\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_r}{r}\frac{\partial a}{\partial r}-\frac{a{\xi }_r{\stackrel{\sim }{\xi }}_r}{r}\frac{\partial \stackrel{\sim }{a}}{\partial r}-{\xi }_r{\stackrel{\sim }{\xi }}_r\frac{\partial a}{\partial r}\frac{\partial \stackrel{\sim }{a}}{\partial r}-\frac{2a\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r}{r}\frac{\partial {\xi }_r}{\partial r}\hfill \\ \hfill & -a\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r\frac{\partial {}^2{\xi }_r}{\partial r^2}-2\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r\frac{\partial {\xi }_r}{\partial r}\frac{\partial a}{\partial r}-\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_r\frac{\partial {}^{2}a}{\partial r^2})S_1\hfill \\ \hfill & +(\frac{a\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r{\xi }_h}{r^2}+\frac{\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r{\xi }_h}{r}\frac{\partial a}{\partial r})(S_2-S_5)-(\frac{\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r{\xi }_h}{r}\frac{\partial a}{\partial r}+\frac{a\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r}{r}\frac{\partial {\xi }_h}{\partial r})\hfill \\ \hfill & ·(S_3+S_6)+(\frac{a\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r{\xi }_h}{r^2}-\frac{a\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r{\xi }_r}{r^2}+\frac{\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r{\xi }_h}{r}\frac{\partial a}{\partial r})S_4\hfill \\ \hfill & +(\frac{a\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_h}{r}\frac{\partial ln{\rho }_0}{\partial r}-\frac{a{\xi }_r{\stackrel{\sim }{\xi }}_h}{r}\frac{\partial \stackrel{\sim }{a}}{\partial r}-\frac{a\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_h}{r}\frac{\partial {\xi }_r}{\partial r}-\frac{\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_h}{r}\frac{\partial a}{\partial r})\hfill \\ \hfill & ·(S_7+S_8)+\frac{a\stackrel{\sim }{a}{\xi }_h{\stackrel{\sim }{\xi }}_h}{r^2}(S_9-S_{10}+S_{11}-2S_{13}+S_{14}-S_{15}\hfill \\ \hfill & -S_{16}-S_{18}-S_{19}-S_{20}-S_{22}-S_{23})-(\frac{a\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_h}{r^2}+\frac{a\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_h}{r}\hfill \\ \hfill & ·\frac{\partial {\xi }_r}{\partial r}+\frac{\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_h}{r}\frac{\partial a}{\partial r})S_{17}-(\frac{a\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_h}{r^2}+\frac{\stackrel{\sim }{a}{\xi }_r{\stackrel{\sim }{\xi }}_h}{r}\frac{\partial a}{\partial r})S_{21}]\text{d}V.\hfill \end{array}$$(50)

Here, we separated the radial contribution from the angular contribution corresponding to the angular kernels S₁ – S₂₃, as given by Equations (G.1) – (G.23). This allowed us to perform the radial and angular integrations independently. We computed the integral over the radial kernels numerically using the composite trapezoidal rule (e.g. Press et al. 2007). For the integrals over the angular kernels, we used the analytical expressions derived in the next section. The radial contributions associated with the angular kernel S₁₂ cancel out. For completeness, however, this kernel is still provided in Appendix G. The general matrix element of a superposition of toroidal magnetic fields, as defined in (45), can be determined by reassigning the indices s and s′ to the two magnetic fields in Equations (46) and (47), respectively, and summing the components.

The linear operator ℱ₀ in the non-magnetic case, given in Equation (7), is Hermitian (Chandrasekhar 1964; Unno et al. 1989). Furthermore, the general force operator in ideal magnetohydrodynamics is self-adjoint and leads to a Hermitian formulation for the total system (Bernstein et al. 1958; Lynden-Bell & Ostriker 1967; Goedbloed & Poedts 2004). However, in our formulation, we decomposed the total force into a modified equilibrium operator and the operator of the direct magnetic effect, as defined in Equation (16): $\mathcal{F}(\mathit{\xi })={\stackrel{\sim }{\mathcal{F}}}_0(\mathit{\xi })+{\mathcal{F}}_B(\mathit{\xi })$. Although the total sum is self-adjoint, the isolated operator of the direct magnetic effect ℱ_B, presented in Equation (17), is not. The terms required to restore the Hermitian symmetry of ℱ_B are located within the pressure terms of the modified equilibrium operator ${\stackrel{\sim }{\mathcal{F}}}_0$. In Appendix C, we demonstrate that our full force operator reduces to the standard self-adjoint expression given in Goedbloed & Poedts (2004) and we explicitly identify the symmetric partners of the ℱ_B terms.

Consequently, the general matrix elements derived solely from ℱ_B do not exhibit Hermitian symmetry (see also Kiefer et al. 2017). More precisely, the matrix is real-valued but asymmetric. The lack of symmetry is structurally evident in terms such as $-a\stackrel{\sim }{a}{\stackrel{\sim }{\xi }}_r\frac{\partial {}^2{\xi }_r}{\partial r^2}$ multiplied by the angular kernel S₁. This term lacks a mirroring pair within the integrand of H_k′k. We verified this numerically by exchanging the modes k′ and k in Equation (50), finding that H_k′k ≠ H_kk′. This asymmetry in principle allows complex eigenvalues and thus complex frequency perturbations. Any complex component of the frequency can be associated with a damping factor for the mode in question. This behaviour contrasts with toroidal flows, such as rotation (Lavely & Ritzwoller 1992), or poloidal flows, such as meridional circulation (Roth & Stix 2008), for which the supermatrix is Hermitian (Kiefer et al. 2017).

4.1. Calculation of integrals over products of scalar generalised spherical harmonics

The angular kernels S₁ – S₂₃ include trigonometric functions and derivatives of scalar GSHs. We represented these derivatives and trigonometric functions in terms of GSHs. We achieved this with the tools derived in Appendix D. The derivatives of GSHs are expressed using Equations (D.5)–(D.9). We eliminated the trigonometric functions using the Equations (D.11), (D.12), (D.14), and (D.15). Equations (G.24)–(G.46) provide the results for the extended angular kernels. Calculating the general matrix element in Equation (50) requires integrals over the angular part of the product of four, five, and six GSHs. We derived complete analytical solutions for these integrals. Previous studies by Dahlen & Tromp (1998) and Kiefer et al. (2017) provide incomplete solutions for these integrals, omitting non-zero contributions.

Dahlen & Tromp (1998) presented a formula for analytically calculating the integrals of the GSHs products on the unit sphere, using the Wigner 3-j symbols described in Appendix E:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{\mathrm{\Omega }}}{Y}_l^{N,m\ast }{Y}_{l_1}^{N_1,m_1}{Y}_{l_2}^{N_2,m_2}\text{d}\mathrm{\Omega }=\&4\pi {(-1)}^{N+m}{\gamma }_l{\gamma }_{l_1}{\gamma }_{l_2}\left(\begin{array}{ccc}l& l_1& l_2\\ -N& N_1& N_2\end{array}\right)\\ \hfill & ·\left(\begin{array}{ccc}l& l_1& l_2\\ -m& m_1& m_2\end{array}\right),\hfill \end{array}$$(51)

where ${\gamma }_l=\sqrt{\frac{2l+1}{4\pi }}$ and the following orthonormality relation is satisfied for GSHs with the same upper index N:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{\mathrm{\Omega }}}{Y}_{l\prime }^{N,m\prime \ast }{Y}_l^{N,m}\text{d}\mathrm{\Omega }={\delta }_{m\prime m}{\delta }_{l\prime l}.\end{array}$$(52)

Applying the orthonormality relation (52) yields the following equation for the product of two scalar GSHs:

$$\begin{array}{cc}\hfill Y_{l_1}^{N_1,m_1}{Y}_{l_2}^{N_2,m_2}=& {\displaystyle \sum _{l,N,m}}4\pi {(-1)}^{N+m}{\gamma }_l{\gamma }_{l_1}{\gamma }_{l_2}\left(\begin{array}{ccc}l& l_1& l_2\\ -N& N_1& N_2\end{array}\right)\hfill \\ & ·\left(\begin{array}{ccc}l& l_1& l_2\\ -m& m_1& m_2\end{array}\right)Y_l^{N,m}.\hfill \end{array}$$(53)

We used the symmetry relation (D.16) to derive a similar result for the product of two generalised Legendre functions:

$$\begin{array}{cc}\hfill P_{l_1}^{N_1,m_1}{P}_{l_2}^{N_2,m_2}=& {\displaystyle \sum _{l,N,m}}(2l+1)\left(\begin{array}{ccc}l& l_1& l_2\\ N& N_1& N_2\end{array}\right)\hfill \\ \hfill & ·\left(\begin{array}{ccc}l& l_1& l_2\\ m& m_1& m_2\end{array}\right)P_l^{N,m}.\hfill \end{array}$$(54)

Equation (51) is valid only when N = N₁ + N₂, because the orthonormality relation (52) applies solely to GSHs with the same upper index N. Equations (51) and (53) facilitate the derivation of analytical solutions for integrals involving products of four or more scalar GSHs. Invoking the selection rule (E.8) for the Wigner 3-j symbols gives the following condition for non-vanishing integrals over products of three or more scalar GSHs:

$$\begin{array}{cc}\hfill m_1=& {\displaystyle \sum _{i>1}}{m}_i.\hfill \end{array}$$(55)

The selection rule (E.8) implies that integrals are non-zero only when N₁ = ∑_i > 1N_i. In the tensorial context considered by Dahlen & Tromp (1998), in which the upper index was associated with the polyads of a tensor field of the appropriate rank, this condition naturally holds, since only then does the integral over the product of GSHs yield a scalar matrix. However, Dahlen & Tromp (1998) imposed that this condition always held, assuming that surface integrals involving the product of two scalar GSHs with N ≠ N′ “never arise in practice”. Our calculations include such integrals over scalar GSHs, which requires the development of a more general method. First, we demonstrated that non-vanishing integrals over the product of three scalar GSHs exist for cases where N ≠ N₁ + N₂. To illustrate this, we evaluated an integral over products of GSHs with low indices. Using $Y_1^{0,1}=-\sqrt{\frac{3}{8\pi }}sin\theta {\text{e}}^{i\varphi }$ and $Y_2^{2,0}=-\sqrt{\frac{5}{4\pi }}\frac{1}{2}\frac{3}{2}{sin}^2\theta$, we find that

$$\begin{array}{c}\hfill {\displaystyle {\int }_{\mathrm{\Omega }}}{Y}_1^{0,1\ast }{Y}_1^{0,1}{Y}_2^{2,0}\text{d}\mathrm{\Omega }=\frac{3}{8}\sqrt{\frac{15}{8\pi }}\frac{16}{15}\ne 0.\end{array}$$

This integral is non-zero while violating the condition N = N₁ + N₂. Moreover, a more general argument follows from considering the symmetry relation (D.17) for even N ≠ 0 and m = 0. In this case, the integral over the product of three scalar GSHs with N = N₁ + N₂ is non-zero according to Equation (51). Using the symmetry relation (D.17) for even N ≠ 0 gives Y_l^N, 0 = Y_l^−N, 0. The integral remains non-zero, but the condition N = N₁ + N₂ is violated. Thus, we argue that assuming integrals over products of scalar GSHs with N₁ ≠ ∑_i > 1N_i vanish, as suggested by Kiefer et al. (2017), is incorrect.

We derived a general solution for integrals over the product of three scalar GSHs. These integrals can be expressed in terms of an integral involving a generalised Legendre function as

$$\begin{array}{cc}& {\displaystyle {\int }_{\mathrm{\Omega }}}{Y}_{l_1}^{N_1,m_1\ast }{Y}_{l_2}^{N_2,m_2}{Y}_{l_3}^{N_3,m_3}\text{d}\mathrm{\Omega }\hfill \\ \hfill & =2\pi {\delta }_{m_1,m_2+m_3}{\gamma }_{l_1}{\gamma }_{l_2}{\gamma }_{l_3}{\displaystyle \sum _{l=|l_2-l_3|}^{l_2+l_3}}{\displaystyle \sum _{l\prime =|l_1-l|}^{l_1+l}}{(4\pi )}^2{\gamma }_l^2{\gamma }_{l\prime }^2\hfill \\ \hfill & ·\left(\begin{array}{ccc}l& l_2& l_3\\ -N_2-N_3& N_2& N_3\end{array}\right)\left(\begin{array}{ccc}l& l_2& l_3\\ -m_1& m_2& m_3\end{array}\right)\hfill \\ \hfill & ·\left(\begin{array}{ccc}l\prime & l_1& l\\ -N_1+N_2+N_3& N_1& -N_2-N_3\end{array}\right)\left(\begin{array}{ccc}l\prime & l_1& l\\ 0& m_1& -m_1\end{array}\right)\hfill \\ \hfill & ·{\displaystyle {\int }_{-1}^1}{P}_{l\prime }^{-N_1+N_2+N_3,0}(\mu )\text{d}\mu ,\hfill \end{array}$$(56)

where Equations (D.1) and (54) are applied, along with the selection rules of the Wigner 3-j symbols (E.8) and (E.10). A closed formula can be found for the remaining integral in Equation (56). With the symmetry relation (D.18), Equation (D.19), and the notation N′= − N₁ + N₂ + N₃, the integral is written as

$$\begin{array}{cc}\hfill {\displaystyle {\int }_{-1}^1}{P}_{l\prime }^{N\prime ,0}(\mu )\text{d}\mu =& {(-1)}^{N\prime }{\displaystyle {\int }_{-1}^1}{P}_{l\prime }^{0,N\prime }(\mu )\text{d}\mu \hfill \\ \hfill =& {(-1)}^{N\prime }\sqrt{\frac{(l-N\prime )!}{(l+N\prime )!}}{\displaystyle {\int }_{-1}^1}{P}_{l\prime }^{N\prime }(\mu )\text{d}\mu .\hfill \end{array}$$(57)

Dong & Lemus (2002) provide a closed formula of the integral over the associated Legendre function P_l^N for N ≥ 0 and l > 0:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{-1}^1}{P}_l^N(\mu )\text{d}\mu =\frac{[{(-1)}^N+{(-1)}^l]2^{N-2}N\mathrm{\Gamma }(\frac{l}{2})\mathrm{\Gamma }(\frac{l+N+1}{2})}{(\frac{l-N}{2})!\mathrm{\Gamma }(\frac{l+3}{2})},\end{array}$$(58)

where Γ(x) denotes the gamma function. The integral is zero if l + N is odd, since then the integrand P_l^N(μ) is an odd function. The symmetry relation (D.21) is applied for the case N < 0. The special case l = 0 is treated separately, because P₀⁰(μ) = 1 simplifies the integral.

Therefore, Equation (56) is fully determined for all |N_i|≤l_i, |m_i|≤l_i, and l_i ≥ 0, where i = 1, 2, 3. Integrals involving the product of four or more GSHs can be reduced to those over three GSHs by repeatedly applying Equation (53). The selection rule (55) is also valid in the general case.

4.2. Selection rules

Two modes, k = (n, l, m) and k′=(n′,l′,m′), couple through the magnetic field if the corresponding general matrix element, which measures the coupling strength between the modes, is non-zero. The general matrix element becomes zero if the integrals over all angular kernels vanish, as indicated by Equation (50). We set the azimuthal orders m of the toroidal magnetic fields to zero. The methods described in Appendix D, which express trigonometric functions in terms of GSHs and compute the derivatives of GSHs, affect only the number of GSHs and their upper index N, leaving the order m unchanged. Hence, Equation (55) is satisfied for all angular kernels if and only if

$$\begin{array}{c}\hfill m=m^{\prime }.\end{array}$$(59)

This result implies the absence of coupling between modes within a multiplet, except for self-coupling.

A second, previously unreported selection rule for mode coupling arises from the symmetries of the angular kernels. All angular kernels exhibit the following symmetry about the equatorial axes:

$$\begin{array}{c}\hfill S_x(\mu )={(-1)}^{l+l\prime +s+s\prime }{S}_x(-\mu ).\end{array}$$(60)

For even l + l′+s + s′, the angular kernels are axially symmetric about the equatorial axis, whereas for odd l + l′+s + s′, the kernels are point symmetric about (θ, S_x(θ)) = (π/2, 0). The integral of the angular kernels with respect to μ = cos θ is symmetric. Consequently, it vanishes if the angular kernels are odd functions of μ. This symmetry yields a zero general matrix element and a zero supermatrix element for odd l + l′+s + s′. Therefore, in addition to Equation (59), modes couple through the toroidal magnetic field only when

$$\begin{array}{c}\hfill l+l\prime +s+s\prime \text{is even}.\end{array}$$(61)

This is especially true for all self-coupling terms in which s = s′.

4.3. Magnetic field model

The radial profiles of the toroidal magnetic field components given by Equations (43) and (44), were modelled using Gaussian functions:

$$\begin{array}{c}\hfill a(r)=\frac{B_{\mathrm{scale}}}{\sigma \sqrt{2\pi }}exp(-0.5(\frac{r-\mu }{\sigma }{)}^2),\end{array}$$(62)

where B_scale scales the distribution to a maximum value B_max, σ represents the Gaussian width, and μ gives the location of the maximum magnetic field strength. Following Antia et al. (2000) derived an upper limit of B_max = 300 kG, we modelled the toroidal magnetic field near the base of the convection zone using B_max = 300 kG, μ = 0.713 R_⊙, σ = 0.05 R_⊙, and harmonic degree s = 2. Figure 1 shows the magnetic field model. The distribution is symmetric in ϕ because the azimuthal orders m of the magnetic field components are set to zero.

Modes are affected by the direct effect of the toroidal magnetic field only if they pass through the layers where the field is present. The modes must penetrate at least to r_thresh = 0.863 R_⊙ to be within 3σ of the magnetic field distribution. The P-modes are confined between a lower turning point r_t and the surface. The lower turning point r_t for p-modes is located where the frequency is equal to the Lamb frequency, ω = S_l(r_t). This gives

$$\begin{array}{c}\hfill \frac{{\omega }^2}{l(l+1)}=\frac{c^2(r_t)}{r_t^2}.\end{array}$$(63)

Multiplets with lower turning points larger than r_thresh are expected to exhibit negligible frequency shifts. Therefore, we excluded multiplets with r_t > r_thresh from this study. This is particularly true for multiplets with l ≥ 138.

5. Results

We computed the frequency shifts in a solar model caused by the direct effect of a toroidal magnetic field located at the base of the convection zone for modes up to 5200 μHz in the range 2 ≤ l ≤ 137. We excluded modes with l ≥ 138, since their lower turning points lie outside the region where the toroidal magnetic field is significant. For each degenerate multiplet ω_0, x² we constructed an eigenspace K_x around ω_ref² = ω_0, x². Contributions from multiplets with frequency gaps larger than about 0.1 μHz are negligible compared to the self-coupling terms. Based on this, we chose an eigenspace width of $\mathrm{\Delta }{\nu }^2=\frac{\mathrm{\Delta }{\omega }^2}{4{\pi }^2}=700\mathrm{\mu }$Hz².

Figure 2 shows the m-averaged frequency shifts as a function of the unperturbed multiplet frequencies. Multiplets with the same harmonic degree l are connected. The overall form of the shifts matches those obtained by Kiefer & Roth (2018) (cf. Figure 3 therein), but their magnitudes, reaching up to 14 nHz, are about twice as large in our study. All multiplets show a positive frequency shift.

Figure 3 presents the m-averaged frequency shifts as a function of the lower turning point r_t. The vertical line marks the radius r_thresh = 0.863 R_⊙, beyond which modes are expected to be unaffected by the direct effect of the toroidal magnetic field. As r_t approaches r_thresh, the frequency shifts of the modes tend to zero. Modes with r_t above r_thresh propagate outside of the magnetised regions, where they experience negligible shifts due to the direct effect of the magnetic field. This can be confirmed by examining the radial contributions in Equation (50). Multiplets with turning points near the maximum of the magnetic field distribution exhibit the greatest frequency shift. The shifts decrease for multiplets with lower turning points beyond the location of the maximum field strength.

Figure 4 displays the lifting of the degeneracy of the multiplet frequency with respect to the azimuthal order m. The absolute frequency shift of the multiplet modes is plotted against m/l and the harmonic degree l. The colour code indicates the radial order n. The degeneracy in m of the unperturbed multiplet is not completely lifted by the interaction with the toroidal magnetic field. Modes within a multiplet remain degenerate if they have the same absolute value of m. This is a result of the symmetry in the angular kernels under the exchange of m ↔ −m, assuming selection rule (59) is fulfilled. The symmetry of the angular kernels is verifiable using Equations (D.1), (D.5), and (D.16). For all multiplets with l ≥ 53 and n ≥ 9, the frequency shift reaches its maximum at m = 0 and decreases monotonically on both sides, approaching zero at |m|=l. For harmonic degrees l ≥ 65, the frequency shift at a fixed m/l increases with radial order n. The mode m = 0, l = 69, n = 19 exhibits the maximum frequency shift within the considered mode range. The plateau of the largest frequency shifts for m = 0 modes aligns with multiplets having lower turning points near the maximum of the magnetic field distribution at μ = 0.713 R_⊙. For multiplets with l ≳ 80, the frequency shift magnitude decreases sharply due to their higher lower turning points. For multiplets with l < 53, only those with intermediate radial order n display a maximum shift at m = 0. The number of multiplets with l < 53 that exhibit this behaviour decreases with decreasing l. For modes in other multiplets, the frequency shifts exhibit a local minimum at m = 0 and two maxima.

6. Discussion and conclusion

Our calculations refine the results of Kiefer et al. (2017) and Kiefer & Roth (2018). While our corrected analysis produces quantitative differences, the qualitative trends and conclusions largely align with their findings. The theory presented here is also applicable to stellar models.

A toroidal magnetic field in the stellar interior couples stellar oscillation eigenmodes and partially lifts the degeneracy of the multiplet frequencies. The general matrix element between two modes quantifies their coupling strength. In this study, we present a revised expression of the general matrix element for the direct effect due to a toroidal magnetic field. Equation (50) is used for forward calculations of the direct effect of a strong toroidal field located at the base of the convection zone. The presented theory enables an investigation of the direct effect of subsurface toroidal magnetic field distributions on stellar oscillation frequencies. We obtained frequency shifts to the oscillation modes by computing the eigenvalues and eigenvectors of the supermatrix following the quasi-degenerate perturbation theory (Lavely & Ritzwoller 1992). Only modes with the same azimuthal order can couple through the toroidal magnetic field (Equation (59)). In addition, the sum of the harmonic degrees of the coupling modes and the toroidal magnetic field components must be even (Equation (61)).

We found that previous studies omit non-zero contributions when calculating integrals over products of scalar GSHs. To address this, we derived a complete method for computing such integrals. The general solution for these integrals, which also includes cases where N₁ ≠ ∑_i > 1N_i, is applicable outside helioseismology, for example in the addition of angular momentum in quantum dynamics.

This study does not include the indirect effect, which represents changes in the equilibrium structure induced by the magnetic field. For magnetic models with s = 2, the indirect effect on the mean multiplet shift is much smaller than the direct effect and can be neglected (Kiefer & Roth 2018). However, a complete description of the frequency perturbations requires its inclusion. Detailed computations of the indirect effect are provided by Kiefer & Roth (2018). The reported absolute values of the indirect effect are likely underestimated because Kiefer & Roth (2018) assumes that integrals with N₁ ≠ ∑_i > 1N_i, which are involved in the computation, vanish. This assumption leads to the neglect of non-zero contributions.

For a strong magnetic field distribution located at the base of the solar convection zone, the m-averaged shift of the multiplet frequencies is of the order of tens of nanohertz. Although this shift is nearly twice that obtained for the same magnetic field model by Kiefer & Roth (2018), it still does not explain the observed frequency shifts over a solar cycle. Broomhall (2017) measured the shift between the maximum and the minimum of solar cycle 23 using data from the Global Oscillation Network Group (GONG). These shifts are of the order of several hundred nanohertz up to microhertz. Hence, near-surface magnetic fields are likely responsible for the observed frequency shifts. However, the existence of deep magnetic field distributions cannot be ruled out. The detection of such fields is strongly hampered, as effects of near-surface fields are stronger by orders of magnitude. The case of the Sun thus differs from that of evolved stars, where a deep fossil field in the stellar core may strongly affect stellar oscillations (Fuller et al. 2015; Bugnet et al. 2021).

The recent study by Vasil et al. (2024), based on magnetohydrodynamic simulations, contributes to the discussion of whether the solar dynamo operates solely near the solar surface. Further studies on the effects of near-surface magnetic fields on oscillations are therefore of particular interest. For this purpose, the solar model should be extended to include an atmosphere. The refined theory is expected to yield significantly different constraints on the model parameters of a near-surface magnetic field that reproduce the observed shifts of Broomhall (2017), compared to those provided by Kiefer & Roth (2018).

Kiefer & Roth (2018) conclude that plotting observed frequency shifts as a function of their lower turning point, as shown in Figure 3, provides hints about the location of the magnetic field configuration. The location of a local maximum in the magnetic field strength can be inferred from the lower turning points of multiplets, whose shifts begin to decline monotonically at a specific harmonic degree. However, the superposition of differently located magnetic field components can hamper such conclusions, since deeper fields produce smaller shifts even if they are stronger than near-surface fields. The revised calculations, which yield larger frequency shifts, combined with sufficient observational resolution, could facilitate identifying dips in such plots that indicate local field maxima.

This work provides a refined analytical solution for the direct effect of a toroidal magnetic field, potentially serving as a stepping stone towards inversion techniques for solar and stellar magnetic field studies. The solution facilitates more detailed studies and tighter constraints on the location of solar and stellar dynamos. Furthermore, we developed a computational framework implementing the quasi-degenerate perturbation theory for a toroidal magnetic field, enabling precise calculations of supermatrices and frequency shifts. This extends its applicability to a wide range of astrophysical scenarios, including different stellar models and magnetic field configurations. For example, it allows parametric studies of detectability thresholds of magnetic field configurations in the Sun and other stars above the measurement noises of related instruments for helio- and asteroseismology (e.g. GONG, HMI, Kepler, TESS, and PLATO). The code is publicly available on GitHub¹, promoting transparency and encouraging community collaboration in future research.

Acknowledgments

The authors thank René Kiefer for fruitful discussions. We thank the anonymous referee for valuable comments that improved the manuscript.

References

---

Appendix A: Mathematical auxiliaries

This appendix lists the vector identities and formulae used to calculate the general matrix element for a toroidal magnetic field described by equation (50). The calculations use spherical coordinates with radius r, co-latitude θ, and azimuthal angle ϕ.

Here, A(r, θ, ϕ) and B(r, θ, ϕ) denote general vector fields. The curl of the cross-product of two vector fields A and B is given by:

$$\begin{array}{c}\hfill \mathrm{\nabla }\times (\mathbf{A}\times \mathbf{B})=\mathbf{A}(\mathrm{\nabla }·\mathbf{B})-\mathbf{B}(\mathrm{\nabla }·\mathbf{A})+(\mathbf{B}·\mathrm{\nabla })\mathbf{A}-(\mathbf{A}·\mathrm{\nabla })\mathbf{B}.\end{array}$$(A.1)

The divergence of a vector field A in spherical coordinates is given by:

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{A}=\frac{1}{r^2}\frac{\partial }{\partial r}(r^2{A}_r)+\frac{1}{rsin\theta }\frac{\partial }{\partial \theta }(A_{\theta }sin\theta )+\frac{1}{rsin\theta }\frac{\partial }{\partial \varphi }{A}_{\varphi },\end{array}$$(A.2)

and its curl is:

$$\begin{array}{cc}\hfill \mathrm{\nabla }\times \mathbf{A}=& \frac{1}{rsin\theta }(\frac{\partial }{\partial \theta }(A_{\varphi }sin\theta )-\frac{\partial A_{\theta }}{\partial \varphi }){\widehat{\mathbf{e}}}_r+\frac{1}{r}(\frac{1}{sin\theta }\frac{\partial A_r}{\partial \varphi }\hfill \\ \hfill & -\frac{\partial (r{A}_{\varphi })}{\partial r}){\widehat{\mathbf{e}}}_{\theta }+\frac{1}{r}(\frac{\partial (r{A}_{\theta })}{\partial r}-\frac{\partial A_r}{\partial \theta }){\widehat{\mathbf{e}}}_{\varphi }.\hfill \end{array}$$(A.3)

The gradient of a scalar field f(r, θ, ϕ) in spherical coordinates is expressed as:

$$\begin{array}{c}\hfill \mathrm{\nabla }f=\frac{\partial f}{\partial r}{\widehat{\mathbf{e}}}_r+\frac{1}{r}\frac{\partial f}{\partial \theta }{\widehat{\mathbf{e}}}_{\theta }+\frac{1}{rsin\theta }\frac{\partial f}{\partial \varphi }{\widehat{\mathbf{e}}}_{\varphi }.\end{array}$$(A.4)

The Laplace operator of a scalar field f in spherical coordinates is defined by:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^{2}f=\frac{1}{r^2}\frac{\partial }{\partial r}(r^2\frac{\partial f}{\partial r})+\frac{1}{r^{2}sin\theta }\frac{\partial }{\partial \theta }(sin\theta \frac{\partial f}{\partial \theta })+\frac{1}{r^2{sin}^2\theta }\frac{\partial {}^{2}f}{\partial {\varphi }^2}.\end{array}$$(A.5)

The material derivative of the vector fields A and B in spherical coordinates can be computed as:

$$\begin{array}{cc}\hfill (\mathbf{A}·\mathrm{\nabla })\mathbf{B}=& (A_r\frac{\partial B_r}{\partial r}+\frac{A_{\theta }}{r}\frac{\partial B_r}{\partial \theta }+\frac{A_{\varphi }}{rsin\theta }\frac{\partial B_r}{\partial \varphi }-\frac{A_{\theta }{B}_{\theta }+A_{\varphi }{B}_{\varphi }}{r}){\widehat{\mathbf{e}}}_r\hfill \\ \hfill & +(A_r\frac{\partial B_{\theta }}{\partial r}+\frac{A_{\theta }}{r}\frac{\partial B_{\theta }}{\partial \theta }+\frac{A_{\varphi }}{rsin\theta }\frac{\partial B_{\theta }}{\partial \varphi }+\frac{A_{\theta }{B}_r}{r}\hfill \\ \hfill & -\frac{A_{\varphi }{B}_{\varphi }cot\theta }{r}){\widehat{\mathbf{e}}}_{\theta }+(A_r\frac{\partial B_{\varphi }}{\partial r}+\frac{A_{\theta }}{r}\frac{\partial B_{\varphi }}{\partial \theta }+\frac{A_{\varphi }}{rsin\theta }\frac{\partial B_{\varphi }}{\partial \varphi }\hfill \\ \hfill & +\frac{A_{\varphi }{B}_r}{r}+\frac{A_{\varphi }{B}_{\theta }cot\theta }{r}){\widehat{\mathbf{e}}}_{\varphi }.\hfill \end{array}$$(A.6)

Appendix B: Perturbation operator of direct effect and calculation of general matrix element

Section 4 presents the general matrix element for the direct effect of a toroidal magnetic field. This appendix outlines the derivation of the general matrix element and provides an explicit expression for the perturbation operator.

Using the mathematical tools described in Appendix A, the operator for a toroidal magnetic field perturbation ℱ_B(ξ_0, k), is derived for fixed values of s and s′ through a straightforward but lengthy calculation. For clarity, the eigenfunctions are represented as follows:

$$\begin{array}{c}\hfill {\mathit{\xi }}_{0,k}=\left(\begin{array}{c}{\xi }^r(r,\theta ,\varphi )\hfill \\ {\xi }^{\theta }(r,\theta ,\varphi )\hfill \\ {\xi }^{\varphi }(r,\theta ,\varphi )\hfill \end{array}\right)=\left(\begin{array}{c}{\xi }_r(r)Y_l^{0,m}\hfill \\ {\xi }_h(r)\frac{\partial }{\partial \theta }{Y}_l^{0,m}\hfill \\ {\xi }_h(r)\frac{1}{sin\theta }\frac{\partial }{\partial \varphi }{Y}_l^{0,m}\hfill \end{array}\right).\end{array}$$(B.1)

The first expressions for the magnetic field components (46) and (47), along with the representation of the eigenfunctions (B.1), are substituted into the magnetic field perturbation operator (17), leading to the following result:

$$\begin{array}{cc}& {\mathcal{F}}_B({\mathit{\xi }}_{0,k})=-\frac{1}{4\pi {\rho }_0}[\{\frac{({\mathit{\xi }}_{0,k}·\mathrm{\nabla }){\rho }_0}{{\rho }_0}+\mathrm{\nabla }·{\mathit{\xi }}_{0,k}\}(\mathrm{\nabla }\times {\mathbf{B}}_s)\times {\mathbf{B}}_{s^{\prime }}\hfill \\ \hfill & +(\mathrm{\nabla }\times {\mathbf{B}}_s)\times \{\mathrm{\nabla }\times ({\mathit{\xi }}_{0,k}\times {\mathbf{B}}_{s^{\prime }})\}+(\mathrm{\nabla }\times \{\mathrm{\nabla }\times ({\mathit{\xi }}_{0,k}\times {\mathbf{B}}_s)\})\times {\mathbf{B}}_{s^{\prime }}]\hfill \\ \hfill & =\frac{1}{4\pi {\rho }_0}[\{(\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r}\frac{\partial ln{\rho }_0}{\partial r}+{\stackrel{\sim }{B}}_{\varphi }{\xi }^r\frac{\partial ln{\rho }_0}{\partial r}\frac{\partial B_{\varphi }}{\partial r}){\widehat{\mathbf{e}}}_r+(\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }{\xi }^{r}cot\theta }{r}\hfill \\ \hfill & ·\frac{\partial ln{\rho }_0}{\partial r}+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r}\frac{\partial ln{\rho }_0}{\partial r}\frac{\partial B_{\varphi }}{\partial \theta }){\widehat{\mathbf{e}}}_{\theta }\}+\{(\frac{2B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r^2}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }{\xi }^{\theta }cot\theta }{r^2}\hfill \\ \hfill & +\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^r}{\partial r}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r^2}\frac{\partial {\xi }^{\theta }}{\partial \theta }+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r^{2}sin\theta }\frac{\partial {\xi }^{\varphi }}{\partial \varphi }+\frac{2{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r}\frac{\partial B_{\varphi }}{\partial r}\hfill \\ \hfill & +\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^{\theta }cot\theta }{r}\frac{\partial B_{\varphi }}{\partial r}+{\stackrel{\sim }{B}}_{\varphi }\frac{\partial {\xi }^r}{\partial r}\frac{\partial B_{\varphi }}{\partial r}+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^{\theta }}{\partial \theta }\frac{\partial B_{\varphi }}{\partial r}+\frac{{\stackrel{\sim }{B}}_{\varphi }}{rsin\theta }\frac{\partial {\xi }^{\varphi }}{\partial \varphi }\hfill \\ \hfill & ·\frac{\partial B_{\varphi }}{\partial r}){\widehat{\mathbf{e}}}_r+(\frac{2B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }{\xi }^{r}cot\theta }{r^2}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }{\xi }^{\theta }{cot}^2\theta }{r^2}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }cot\theta }{r}\frac{\partial {\xi }^r}{\partial r}\hfill \\ \hfill & +\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }cot\theta }{r^2}\frac{\partial {\xi }^{\theta }}{\partial \theta }+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }cot\theta }{r^{2}sin\theta }\frac{\partial {\xi }^{\varphi }}{\partial \varphi }+\frac{2{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r^2}\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^{\theta }cot\theta }{r^2}\hfill \\ \hfill & ·\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^r}{\partial r}\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r^2}\frac{\partial {\xi }^{\theta }}{\partial \theta }\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r^{2}sin\theta }\frac{\partial {\xi }^{\varphi }}{\partial \varphi }\frac{\partial B_{\varphi }}{\partial \theta }){\widehat{\mathbf{e}}}_{\theta }\}\hfill \\ \hfill & -\{(\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r^2}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^r}{\partial r}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r^2}\frac{\partial {\xi }^{\theta }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r}\frac{\partial B_{\varphi }}{\partial r}+\frac{B_{\varphi }{\xi }^r}{r}\frac{\partial {\stackrel{\sim }{B}}_{\varphi }}{\partial r}\hfill \\ \hfill & +\frac{B_{\varphi }{\xi }^{\theta }}{r^2}\frac{\partial {\stackrel{\sim }{B}}_{\varphi }}{\partial \theta }+{\stackrel{\sim }{B}}_{\varphi }\frac{\partial {\xi }^r}{\partial r}\frac{\partial B_{\varphi }}{\partial r}+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^{\theta }}{\partial \theta }\frac{\partial B_{\varphi }}{\partial r}+{\xi }^r\frac{\partial B_{\varphi }}{\partial r}\frac{\partial {\stackrel{\sim }{B}}_{\varphi }}{\partial r}+\frac{{\xi }^{\theta }}{r}\hfill \\ \hfill & ·\frac{\partial B_{\varphi }}{\partial r}\frac{\partial {\stackrel{\sim }{B}}_{\varphi }}{\partial \theta }){\widehat{\mathbf{e}}}_r+(\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }{\xi }^{r}cot\theta }{r^2}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }cot\theta }{r}\frac{\partial {\xi }^r}{\partial r}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }cot\theta }{r^2}\hfill \\ \hfill & ·\frac{\partial {\xi }^{\theta }}{\partial \theta }+\frac{B_{\varphi }{\xi }^{r}cot\theta }{r}\frac{\partial {\stackrel{\sim }{B}}_{\varphi }}{\partial r}+\frac{B_{\varphi }{\xi }^{\theta }cot\theta }{r^2}\frac{\partial {\stackrel{\sim }{B}}_{\varphi }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r^2}\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^r}{\partial r}\hfill \\ \hfill & ·\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r^2}\frac{\partial {\xi }^{\theta }}{\partial \theta }\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\xi }^r}{r}\frac{\partial {\stackrel{\sim }{B}}_{\varphi }}{\partial r}\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\xi }^{\theta }}{r^2}\frac{\partial {\stackrel{\sim }{B}}_{\varphi }}{\partial \theta }\frac{\partial B_{\varphi }}{\partial \theta }){\widehat{\mathbf{e}}}_{\theta }\hfill \\ \hfill & +(\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r^{2}sin\theta }\frac{\partial {\xi }^r}{\partial \varphi }+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }cot\theta }{r^{2}sin\theta }\frac{\partial {\xi }^{\theta }}{\partial \varphi }+\frac{{\stackrel{\sim }{B}}_{\varphi }}{rsin\theta }\frac{\partial {\xi }^r}{\partial \varphi }\frac{\partial B_{\varphi }}{\partial r}+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r^{2}sin\theta }\hfill \\ \hfill & ·\frac{\partial {\xi }^{\theta }}{\partial \varphi }\frac{\partial B_{\varphi }}{\partial \theta }){\widehat{\mathbf{e}}}_{\varphi }\}-\{(\frac{2B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^r}{\partial r}+\frac{2{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r}\frac{\partial B_{\varphi }}{\partial r}+B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }\frac{\partial {}^2{\xi }^r}{\partial r^2}\hfill \\ \hfill & +\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r^2{sin}^2\theta }\frac{\partial {}^2{\xi }^r}{\partial {\varphi }^2}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {}^2{\xi }^{\theta }}{\partial r\partial \theta }+2{\stackrel{\sim }{B}}_{\varphi }\frac{\partial {\xi }^r}{\partial r}\frac{\partial B_{\varphi }}{\partial r}+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^{\theta }}{\partial \theta }\frac{\partial B_{\varphi }}{\partial r}\hfill \\ \hfill & +\frac{{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^{\theta }}{\partial r}\frac{\partial B_{\varphi }}{\partial \theta }+{\xi }^r{\stackrel{\sim }{B}}_{\varphi }\frac{\partial {}^2{B}_{\varphi }}{\partial r^2}+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^{\theta }}{r}\frac{\partial {}^2{B}_{\varphi }}{\partial r\partial \theta }){\widehat{\mathbf{e}}}_r+(\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }{\xi }^{r}cot\theta }{r^2}\hfill \\ \hfill & +\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r^2}\frac{\partial {\xi }^r}{\partial \theta }+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }cot\theta }{r}\frac{\partial {\xi }^r}{\partial r}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }cot\theta }{r^2}\frac{\partial {\xi }^{\theta }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^{r}cot\theta }{r}\hfill \\ \hfill & ·\frac{\partial B_{\varphi }}{\partial r}+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r^2}\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^{\theta }cot\theta }{r^2}\frac{\partial B_{\varphi }}{\partial \theta }+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r^2}\frac{\partial {}^2{\xi }^{\theta }}{\partial {\theta }^2}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r^2{sin}^2\theta }\hfill \\ \hfill & ·\frac{\partial {}^2{\xi }^{\theta }}{\partial {\varphi }^2}+\frac{B_{\varphi }{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {}^2{\xi }^r}{\partial r\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^r}{\partial r}\frac{\partial B_{\varphi }}{\partial \theta }+\frac{2{\stackrel{\sim }{B}}_{\varphi }}{r^2}\frac{\partial {\xi }^{\theta }}{\partial \theta }\frac{\partial B_{\varphi }}{\partial \theta }+\frac{{\stackrel{\sim }{B}}_{\varphi }}{r}\frac{\partial {\xi }^r}{\partial \theta }\hfill \\ \hfill & ·\frac{\partial B_{\varphi }}{\partial r}+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^{\theta }}{r^2}\frac{\partial {}^2{B}_{\varphi }}{\partial {\theta }^2}+\frac{{\stackrel{\sim }{B}}_{\varphi }{\xi }^r}{r}\frac{\partial {}^2{B}_{\varphi }}{\partial r\partial \theta }){\widehat{\mathbf{e}}}_{\theta }\}].\hfill \end{array}$$(B.2)

The four terms in the first equality are kept separate and enclosed within curly brackets. The left scalar product of equation (B.2) with ρ₀ξ_0, k′^* is computed and integrated over the entire stellar volume. Subsequently, the second expressions from equations (46), (47), and (B.1) are inserted, using the notation for the complex conjugated eigenfunctions given in (49). Finally, sorting by the angular kernels S₁ – S₂₃, as defined in equations (G.1) – (G.23) in Appendix G, results in the general matrix element for the direct effect of a toroidal magnetic field, given in equation (50).

Appendix C: Restoration of hermiticity in the total force operator

In this appendix, we demonstrate that our full force operator, defined in equation (16), reduces to the expression given in equation (6.23) of Goedbloed & Poedts (2004). We note that Goedbloed & Poedts (2004) used a magnetised equilibrium, did not include perturbations of the gravitational potential, and did not use cgs units.

We begin by rewriting the first term of the operator of the direct magnetic effect ℱ_B, presented in equation (17), as:

$$\begin{array}{cc}& -\frac{1}{4\pi {\stackrel{\sim }{\rho }}_0}\{\frac{(\mathit{\xi }·\mathrm{\nabla }){\stackrel{\sim }{\rho }}_0}{{\stackrel{\sim }{\rho }}_0}+\mathrm{\nabla }·\mathit{\xi }\}(\mathrm{\nabla }\times {\stackrel{\sim }{\mathbf{B}}}_0)\times {\stackrel{\sim }{\mathbf{B}}}_0\hfill \\ \hfill =& -\{\frac{(\mathit{\xi }·\mathrm{\nabla }){\stackrel{\sim }{\rho }}_0+{\stackrel{\sim }{\rho }}_0\mathrm{\nabla }·\mathit{\xi }}{{\stackrel{\sim }{\rho }}_0^2}\}(\mathrm{\nabla }{\stackrel{\sim }{p}}_0+{\stackrel{\sim }{\rho }}_0\mathrm{\nabla }{\stackrel{\sim }{\varphi }}_0)\hfill \\ \hfill =& -\frac{1}{{\stackrel{\sim }{\rho }}_0^2}\mathrm{\nabla }·({\stackrel{\sim }{\rho }}_0\mathit{\xi })(\mathrm{\nabla }{\stackrel{\sim }{p}}_0+{\stackrel{\sim }{\rho }}_0\mathrm{\nabla }{\stackrel{\sim }{\varphi }}_0)\hfill \\ \hfill =& \frac{{\rho }^{\prime }}{{\stackrel{\sim }{\rho }}_0^2}\mathrm{\nabla }{\stackrel{\sim }{p}}_0+\frac{{\rho }^{\prime }}{{\stackrel{\sim }{\rho }}_0}\mathrm{\nabla }{\stackrel{\sim }{\varphi }}_0.\hfill \end{array}$$(C.1)

Here, we used the magnetised equilibrium condition (14) to eliminate the term containing the equilibrium magnetic field. We also applied the product rule for divergence:

$$\begin{array}{c}\hfill \mathrm{\nabla }·(f\mathit{A})=\mathit{A}·\mathrm{\nabla }f+f\mathrm{\nabla }·\mathit{A},\end{array}$$(C.2)

and the linearised continuity equation (e.g., Lynden-Bell & Ostriker 1967; Christensen-Dalsgaard 2002; Goedbloed & Poedts 2004):

$$\begin{array}{c}\hfill {\rho }^{\prime }=-\mathrm{\nabla }({\stackrel{\sim }{\rho }}_0\mathit{\xi }).\end{array}$$(C.3)

We can express the modified equilibrium operator ${\stackrel{\sim }{\mathcal{F}}}_0$ as:

$$\begin{array}{cc}\hfill {\stackrel{\sim }{\mathcal{F}}}_0(\mathit{\xi })& =\frac{1}{{\stackrel{\sim }{\rho }}_0}\mathrm{\nabla }{p}^{\prime }-{\mathit{g}}^{\prime }-\frac{{\rho }^{\prime }}{{\stackrel{\sim }{\rho }}_0^2}\mathrm{\nabla }{\stackrel{\sim }{p}}_0\hfill \end{array}$$(C.4)

$$\begin{array}{cc}& =-\frac{1}{{\stackrel{\sim }{\rho }}_0}\mathrm{\nabla }(\mathit{\xi }·\mathrm{\nabla }{\stackrel{\sim }{p}}_0+\gamma {\stackrel{\sim }{p}}_0\mathrm{\nabla }·\mathit{\xi })-{\mathit{g}}^{\prime }+\frac{\mathrm{\nabla }·({\stackrel{\sim }{\rho }}_0\mathit{\xi })}{{\stackrel{\sim }{\rho }}_0^2}\mathrm{\nabla }{\stackrel{\sim }{p}}_0,\hfill \end{array}$$(C.5)

where we used equation (C.3) and the linearised energy equation (e.g., Lynden-Bell & Ostriker 1967; Christensen-Dalsgaard 2002; Goedbloed & Poedts 2004):

$$\begin{array}{c}\hfill p^{\prime }=-\mathit{\xi }·\mathrm{\nabla }{\stackrel{\sim }{p}}_0-\gamma {\stackrel{\sim }{p}}_0\mathrm{\nabla }·\mathit{\xi }.\end{array}$$(C.6)

Plugging equation (C.1) into the expression for ℱ_B given in equation (17) and summing it with ${\stackrel{\sim }{\mathcal{F}}}_0$ yields:

$$\begin{array}{cc}\hfill \mathcal{F}& ={\stackrel{\sim }{\mathcal{F}}}_0+{\mathcal{F}}_B\hfill \\ \hfill & =\frac{1}{{\stackrel{\sim }{\rho }}_0}\mathrm{\nabla }{p}^{\prime }+\frac{{\rho }^{\prime }}{{\stackrel{\sim }{\rho }}_0}\mathrm{\nabla }{\stackrel{\sim }{\varphi }}_0+\mathrm{\nabla }{\varphi }^{\prime }\hfill \\ \hfill & -\frac{1}{4\pi {\stackrel{\sim }{\rho }}_0}((\mathrm{\nabla }\times {\stackrel{\sim }{\mathbf{B}}}_0)\times {\mathbf{B}}^{\prime }+(\mathrm{\nabla }\times {\mathbf{B}}^{\prime })\times {\stackrel{\sim }{\mathbf{B}}}_0).\hfill \end{array}$$(C.7)

This is equivalent to $-\frac{1}{\stackrel{\sim }{\rho }}{\mathbf{F}}_{\mathrm{GP}}(\mathit{\xi })$, where F_GP is given in equation (6.23) of Goedbloed & Poedts (2004), except that we applied the linearised continuity equation and perturbed the gravitational potential resulting in an additional ∇ϕ′ term.

Next, we explicitly identify where the Hermiticity breaks down in ℱ_B and find the symmetric partners in ${\stackrel{\sim }{\mathcal{F}}}_0$. We use the inner product:

$$\begin{array}{c}\hfill ⟨\mathit{\eta },\mathit{\xi }⟩={\displaystyle {\int }_V}{\mathit{\eta }}^{\ast }·\mathit{\xi }{\stackrel{\sim }{\rho }}_0\mathrm{d}V,\end{array}$$(C.8)

where η, ξ are two displacement fields. To simplify the proof and follow the formalism of Goedbloed & Poedts (2004) more closely, we assumed fixed boundary conditions:

$$\begin{array}{c}\hfill \mathit{\xi }·\mathit{n}=0\text{at}r=R,\end{array}$$(C.9)

where n is the normal vector of the surface. Under these boundary conditions, surface integrals vanish; these are equivalent to the volume integrals of divergence terms after applying the Gauss-divergence theorem. For the zero-boundary conditions (Unno et al. 1989) discussed in Section 2.1, the surface terms do not vanish but instead form a symmetric boundary energy term (Bernstein et al. 1958; Lynden-Bell & Ostriker 1967). Therefore, Hermiticity is preserved even when ξ ⋅ n ≠ 0. We split ℱ_B(ξ) into three parts:

$$\begin{array}{cc}\hfill \mathcal{M}(\mathit{\xi })& =-\frac{1}{4\pi {\stackrel{\sim }{\rho }}_0}((\mathrm{\nabla }\times {\stackrel{\sim }{\mathbf{B}}}_0)\times {\mathbf{B}}^{\prime }),\hfill \end{array}$$(C.10)

$$\begin{array}{cc}\hfill \mathcal{N}(\mathit{\xi })& =-\frac{1}{4\pi {\stackrel{\sim }{\rho }}_0}((\mathrm{\nabla }\times {\mathbf{B}}^{\prime })\times {\stackrel{\sim }{\mathbf{B}}}_0),\hfill \end{array}$$(C.11)

$$\begin{array}{cc}\hfill \mathcal{O}(\mathit{\xi })& =\frac{{\rho }^{\prime }}{{\stackrel{\sim }{\rho }}_0^2}\mathrm{\nabla }{\stackrel{\sim }{p}}_0+\frac{{\rho }^{\prime }}{{\stackrel{\sim }{\rho }}_0}\mathrm{\nabla }{\stackrel{\sim }{\varphi }}_0,\hfill \end{array}$$(C.12)

where in the last operator we applied equation (C.1). Goedbloed & Poedts (2004) show that ${\mathit{\eta }}^{\ast }·\mathcal{N}(\mathit{\xi }){\stackrel{\sim }{\rho }}_0$ and the second pressure term ${\mathit{\eta }}^{\mathbf{\ast }}·\mathrm{\nabla }(\gamma {\stackrel{\sim }{p}}_0\mathrm{\nabla }·\mathit{\xi })$ form symmetric expressions plus divergence terms, which can be discarded upon integration due to the boundary condition (C.9). We can also show explicitly that the work done by the perturbed gravitational field is Hermitian:

$$\begin{array}{cc}\hfill -{\displaystyle {\int }_V}{\mathit{\eta }}^{\ast }·{\mathit{g}}^{\prime }(\mathit{\xi }){\stackrel{\sim }{\rho }}_0\mathrm{d}V=& {\displaystyle {\int }_V}{\mathit{\eta }}^{\ast }·\mathrm{\nabla }{\varphi }^{\prime }[\mathit{\xi }]{\stackrel{\sim }{\rho }}_0\mathrm{d}V\hfill \\ \hfill \stackrel{(\mathrm{C}.2)}{=}& {\displaystyle {\int }_V}\mathrm{\nabla }·({\mathit{\eta }}^{\ast }{\varphi }^{\prime }[\mathit{\xi }]{\stackrel{\sim }{\rho }}_0)\mathrm{d}V\hfill \\ \hfill & -{\displaystyle {\int }_V}{\varphi }^{\prime }[\mathit{\xi }]\mathrm{\nabla }·({\stackrel{\sim }{\rho }}_0{\mathit{\eta }}^{\ast })\mathrm{d}V,\hfill \end{array}$$(C.13)

where we used g′(ξ) =  − ∇ϕ′[ξ] in the first equality. The first integral in the last equation is again a vanishing boundary term. Using the continuity equation in the form $\rho \prime [{\mathit{\eta }}^{*}]=-\mathrm{\nabla }·({\stackrel{\sim }{\rho }}_0{\mathit{\eta }}^{*})$ and the Poisson equation ∇²ϕ′ = 4πGρ′, we get:

$$\begin{array}{cc}\hfill -{\displaystyle {\int }_V}{\mathit{\eta }}^{\ast }·{\mathit{g}}^{\prime }(\mathit{\xi }){\rho }_0\mathrm{d}V=& {\displaystyle {\int }_V}{\varphi }^{\prime }[\mathit{\xi }]{\rho }^{\prime }[{\mathit{\eta }}^{\mathbf{\ast }}]\mathrm{d}V\hfill \\ \hfill =& \frac{1}{4\pi G}{\displaystyle {\int }_V}{\varphi }^{\prime }[\mathit{\xi }]{\mathrm{\nabla }}^2{\varphi }^{\prime }[{\mathit{\eta }}^{\mathbf{\ast }}]\mathrm{d}V\hfill \\ \hfill \stackrel{(\mathrm{C}.2)}{=}& \frac{1}{4\pi G}[{\displaystyle {\int }_V}\mathrm{\nabla }·({\varphi }^{\prime }[\mathit{\xi }]\mathrm{\nabla }{\varphi }^{\prime }[{\mathit{\eta }}^{\mathbf{\ast }}])\mathrm{d}V\hfill \\ \hfill & -{\displaystyle {\int }_V}\mathrm{\nabla }{\varphi }^{\prime }[\mathit{\xi }]·\mathrm{\nabla }{\varphi }^{\prime }[{\mathit{\eta }}^{\mathbf{\ast }}]\mathrm{d}V]\hfill \\ \hfill =& -\frac{1}{4\pi G}{\displaystyle {\int }_V}\mathrm{\nabla }{\varphi }^{\prime }[\mathit{\xi }]·\mathrm{\nabla }{\varphi }^{\prime }[{\mathit{\eta }}^{\mathbf{\ast }}]\mathrm{d}V,\hfill \end{array}$$(C.14)

where in the last step we again discard the boundary term.

The remaining terms are exactly the inhomogeneity terms given in Goedbloed & Poedts (2004):

$$\begin{array}{cc}& {\mathit{\eta }}^{\mathbf{\ast }}[-\mathrm{\nabla }(\mathit{\xi }·\mathrm{\nabla }{\stackrel{\sim }{p}}_0)-\frac{{\rho }^{\prime }}{{\stackrel{\sim }{\rho }}_0}\mathrm{\nabla }{\stackrel{\sim }{p}}_0+\mathcal{O}(\mathit{\xi })+\mathcal{M}(\mathit{\xi })]\hfill \\ \hfill =& -{\mathit{\eta }}^{\mathbf{\ast }}[\mathrm{\nabla }(\mathit{\xi }·\mathrm{\nabla }{\stackrel{\sim }{p}}_0)+\mathrm{\nabla }·({\stackrel{\sim }{\rho }}_0\mathit{\xi })\mathrm{\nabla }{\stackrel{\sim }{\varphi }}_0+\frac{1}{4\pi }(\mathrm{\nabla }\times {\stackrel{\sim }{\mathbf{B}}}_0)\times {\mathbf{B}}^{\prime }].\hfill \end{array}$$(C.15)

Following the derivation in Goedbloed & Poedts (2004), these terms reduce to a symmetric form plus divergence terms that vanish upon integration. Thus, we have shown that the Hermiticity of ℳ(ξ) and 𝒪(ξ) is restored by combining them with terms in ${\stackrel{\sim }{\mathcal{F}}}_0$.

Appendix D: Properties of generalised spherical harmonics

Here, an overview of generalised spherical harmonics (GSHs) and their properties is presented based on the conventions from Dahlen & Tromp (1998). The GSHs constitute an orthonormal basis, allowing tensor fields of any order to be expanded in terms of GSHs. Scalar GSHs are used to represent the oscillation eigenfunctions and the toroidal magnetic fields in this study. This approach facilitates the analytical calculation of integrals over products of scalar GSHs. All essential identities and relations of the GSHs and the generalised Legendre functions necessary for our calculations are provided. The formulae from Dahlen & Tromp (1998) are applied to derive expressions for derivatives of GSHs and eliminate trigonometric functions.

The scalar GSHs Y_l^N, m(θ, ϕ) are defined as:

$$\begin{array}{c}\hfill Y_l^{N,m}(\theta ,\varphi )={\gamma }_l{P}_l^{N,m}(\mu ){\text{e}}^{im\varphi },\end{array}$$(D.1)

where μ = cos θ, ${\gamma }_l=\sqrt{\frac{2l+1}{4\pi }}$, and P_l^N, m(μ) are the generalised Legendre functions. In this notation, the GSH reduces to the ordinary spherical harmonic function Y_l^0, m ≡ Y_l^m for N = 0. The generalised Legendre functions of degree l, order −l ≤ m ≤ l, and upper index −l ≤ N ≤ l are defined by Phinney & Burridge (1973) as:

$$\begin{array}{cc}\hfill P_l^{N,m}=& \frac{1}{2^l}\sqrt{\frac{(l+m)!}{(l+N)!(l-N)!(l-m)!}}(1-\mu {)}^{-\frac{1}{2}(m-N)}\hfill \\ \hfill & ·(1+\mu {)}^{-\frac{1}{2}(m+N)}\frac{{\text{d}}^{l-m}}{\text{d}{\mu }^{l-m}}[(\mu -1{)}^{l-N}(\mu +1{)}^{l+N}],\hfill \end{array}$$(D.2)

and satisfy the following ordinary differential equation:

$$\begin{array}{cc}& \{(1-{\mu }^2)\frac{{\text{d}}^2}{\text{d}{\mu }^2}-2\mu \frac{\text{d}}{\text{d}\mu }+[l(l+1)\hfill \\ \hfill & -\frac{m^2-2mN\mu +N^2}{1-{\mu }^2}]\}{P}_l^{N,m}(\mu )=0.\hfill \end{array}$$(D.3)

This differential equation has non-zero solutions in the form of equation (D.2) only for |m|≤l and |N|≤l, while it is trivially solved for other values. Therefore, a logical extension for the definition of the generalised Legendre functions is to set P_l^N, m = 0 for |N|> l or |m|> l.

For a fixed N, the generalised Legendre functions fulfil the following recursion relation:

$$\begin{array}{cc}\hfill (\frac{N}{sin\theta }-mcot\theta )P_l^{N,m}(\mu )=& \frac{1}{2}(\sqrt{(l+m)(l-m+1)}{P}_l^{N,m-1}\hfill \\ \hfill & +\sqrt{(l-m)(l+m+1)}{P}_l^{N,m+1}).\hfill \end{array}$$(D.4)

This equation forms the foundation for numerically computing the (2l + 1)⋅(2l + 1) generalised Legendre functions for a given l, with |m|≤l and |N|≤l (e.g. Masters & Richards-Dinger 1998). The longitudinal derivative of order k of a GSH, as derived from equation (D.1), is:

$$\begin{array}{c}\hfill \frac{\partial {}^k}{\partial {\varphi }^k}{Y}_l^{N,m}={(im)}^k{Y}_l^{N,m}.\end{array}$$(D.5)

The co-latitudinal derivative of a GSH is expressed in terms of Y_l^{N ± 1, m} as:

$$\begin{array}{c}\hfill \frac{\partial }{\partial \theta }{Y}_l^{N,m}=\frac{1}{2}(c_l^N{Y}_l^{N-1,m}-c_l^{N+1}{Y}_l^{N+1,m}),\end{array}$$(D.6)

where $c_l^N=\sqrt{(l+N)(l-N+1)}$ is a constant that satisfies two symmetry properties:

$$\begin{array}{c}\hfill c_l^0=c_l^1\text{and}{c}_l^{-N}=c_l^{N+1}.\end{array}$$(D.7)

We regularly employ these symmetry relations throughout the thesis. To obtain higher-order derivatives, we apply equation (D.6) iteratively. For example, the second- and third-order co-latitudinal derivatives are given by:

$$\begin{array}{cc}\hfill \frac{\partial {}^2}{\partial {\theta }^2}{Y}_l^{N,m}=& \frac{1}{4}[c_l^N{c}_l^{N-1}{Y}_l^{N-2,m}-(c_l^N{c}_l^N+c_l^{N+1}{c}_l^{N+1})Y_l^{N,m}\hfill \\ \hfill & +c_l^{N+1}{c}_l^{N+2}{Y}_l^{N+2,m}],\hfill \end{array}$$(D.8)

$$\begin{array}{cc}\hfill \frac{\partial {}^3}{\partial {\theta }^3}{Y}_l^{N,m}=& \frac{1}{8}[c_l^N{c}_l^{N-1}{c}_l^{N-2}{Y}_l^{N-3,m}-(c_l^N{c}_l^{N-1}{c}_l^{N-1}+c_l^N{c}_l^N{c}_l^N\hfill \\ \hfill & +c_l^N{c}_l^{N+1}{c}_l^{N+1})Y_l^{N-1,m}+(c_l^N{c}_l^N{c}_l^{N+1}+c_l^{N+1}{c}_l^{N+1}{c}_l^{N+1}\hfill \\ \hfill & +c_l^{N+1}{c}_l^{N+2}{c}_l^{N+2})Y_l^{N+1,m}-c_l^{N+1}{c}_l^{N+2}{c}_l^{N+3}{Y}_l^{N+3,m}].\hfill \end{array}$$(D.9)

We discuss the cases in which these recursive relations are valid. In general, we define GSHs only for |N|≤l. Therefore, we conclude that |N|≤l − k, where k is the order of the derivative. However, by setting Y_l^N, m = 0 for |N|> l, the relations are valid for |N|≤l − k + 1, as the prefactors of the terms containing a GSH with |N|=l − k + 1 are zero. Thus, the value assigned to Y_l^N, m for |N|> l is arbitrary. In other cases, the formulae break down because some of the prefactors become imaginary. Therefore, we can also use equation (D.9) for l = 2 and N = 0. This is relevant for the calculation with s = 2 in the angular kernel S₂₀, given in equation (G.20). Similar reasoning applies to the recursion relations discussed later in this section.

The GSHs satisfy the following recursion relation for |N|≤l:

$$\begin{array}{c}\hfill [Ncot\theta -\frac{m}{sin\theta }]Y_l^{N,m}=\frac{1}{2}(c_l^N{Y}_l^{N-1,m}+c_l^{N+1}{Y}_l^{N+1,m}).\end{array}$$(D.10)

We aim to express trigonometric functions in terms of GSHs in order to exploit methods for calculating integrals over products of GSHs. We can find suitable relations by setting N = 0 or m = 0 in the recursion relation above. This yields a relation for eliminating a factor of $\frac{1}{sin\theta }$ when N = 0, m ≠ 0, and l ≥ 1:

$$\begin{array}{c}\hfill \frac{1}{sin\theta }{Y}_l^{0,m}=-\frac{1}{2m}{c}_l^0(Y_l^{-1,m}+Y_l^{1,m}).\end{array}$$(D.11)

If m = 0, N ≠ 0, and |N|≤l, then we obtain an expression for cotθ:

$$\begin{array}{c}\hfill cot\theta Y_l^{N,0}=\frac{1}{2N}(c_l^N{Y}_l^{N-1,0}+c_l^{N+1}{Y}_l^{N+1,0}).\end{array}$$(D.12)

We can derive further relations by taking the co-latitudinal derivative of the recursion relation (D.10):

$$\begin{array}{cc}& [Ncot\theta -\frac{m}{sin\theta }]\frac{1}{2}(c_l^N{Y}_l^{N-1,m}+c_l^{N+1}{Y}_l^{N+1,m})\hfill \\ \hfill & +[\frac{mcot\theta }{sin\theta }-\frac{N}{{sin}^2\theta }]Y_l^{N,m}\hfill \\ \hfill & =\frac{1}{4}[c_l^N{c}_l^{N-1}{Y}_l^{N-2,m}+(c_l^{N+1}{c}_l^{N+1}-c_l^N{c}_l^N)Y_l^{N,m}\hfill \\ \hfill & -c_l^{N+1}{c}_l^{N+2}{Y}_l^{N+2,m}].\hfill \end{array}$$(D.13)

Substituting N = 0 into this equation and using the identity $Y_1^{1,0}={\gamma }_1{P}_1^{1,0}(\mu )=\sqrt{\frac{3}{8\pi }}sin\theta$ gives the following for m ≠ 0 and l ≥ 1:

$$\begin{array}{cc}\hfill cot\theta Y_l^{0,m}=& \frac{1}{2}{c}_l^0(Y_l^{-1,m}-Y_l^{1,m})+\frac{1}{m}\sqrt{\frac{\pi }{6}}{c}_l^0{c}_l^2(Y_l^{-2,m}{Y}_1^{1,0}\hfill \\ \hfill & -Y_l^{2,m}{Y}_1^{1,0}).\hfill \end{array}$$(D.14)

If we set m = 0 and use in addition to $Y_1^{1,0}=\sqrt{\frac{3}{8\pi }}sin\theta$ the identity $Y_1^{0,0}=\sqrt{\frac{3}{4\pi }}cos\theta$, we get for N ≠ 0 and |N|≤l − 1:

$$\begin{array}{cc}\hfill \frac{1}{sin\theta }{Y}_l^{N,0}=& \sqrt{\frac{\pi }{3}}(c_l^N{Y}_l^{N-1,0}{Y}_1^{0,0}-c_l^{N+1}{Y}_l^{N+1,0}{Y}_1^{0,0})\hfill \\ \hfill & -\frac{1}{N}\sqrt{\frac{\pi }{6}}[c_l^N{c}_l^{N-1}{Y}_l^{N-2,0}{Y}_1^{1,0}+(c_l^{N+1}{c}_l^{N+1}\hfill \\ \hfill & -c_l^N{c}_l^N)Y_l^{N,0}{Y}_1^{1,0}-c_l^{N+1}{c}_l^{N+2}{Y}_l^{N+2,0}{Y}_1^{1,0}].\hfill \end{array}$$(D.15)

We emphasise that this relation remains valid even for l = 2 and N = 1. This is particularly important for calculating the angular kernels S₁₁, S₁₄, S₁₉, S₂₂, and S₂₃, given by equations (G.11), (G.14), (G.19), (G.22), and (G.23), respectively, where s = s′ = 2. Therefore, no modifications to the equations are required for the case of s = s′ = 2.

An important symmetry relation for the GSHs is derived from the following symmetry relation of the generalised Legendre function:

$$\begin{array}{c}\hfill P_l^{-N,-m}(\mu )={(-1)}^{m+N}{P}_l^{N,m}(\mu ).\end{array}$$(D.16)

For m = 0, the GSHs satisfy:

$$\begin{array}{c}\hfill Y_l^{-N,0}(\mu ,\varphi )={(-1)}^N{Y}_l^{N,0}(\mu ,\varphi ).\end{array}$$(D.17)

This relation is applied throughout the calculation of the explicit expressions of the angular kernels. Additionally, we have the following symmetry relation for the generalised Legendre functions:

$$\begin{array}{cc}\hfill P_l^{m,N}(\mu )& ={(-1)}^{m+N}{P}_l^{N,m}(\mu ).\hfill \end{array}$$(D.18)

This relation is used to derive the solution to the integral over the product of three GSHs.

We relate the generalised Legendre functions to the associated Legendre functions P_l^m(μ):

$$\begin{array}{c}\hfill P_l^{0,m}(\mu )=\sqrt{\frac{(l-m)!}{(l+m)!}}{P}_l^m(\mu ),\end{array}$$(D.19)

where we define the associated Legendre function in contrast to Dahlen & Tromp (1998) as (Arfken et al. 2012):

$$\begin{array}{c}\hfill P_l^m(\mu )=\frac{{(-1)}^m}{2^{l}l!}(1-{\mu }^2{)}^{m/2}\frac{{\text{d}}^{l+m}}{\text{d}{\mu }^{l+m}}({\mu }^2-1{)}^l.\end{array}$$(D.20)

This definition is more commonly used. This is important for the correct application of the symmetric integral formula over the associated Legendre function, as given in equation (58).

The associated Legendre functions obey the following symmetry:

$$\begin{array}{c}\hfill P_l^{-m}={(-1)}^m\frac{(l-m)!}{(l+m)!}{P}_l^m.\end{array}$$(D.21)

The following three GSH relations are used to analyse the symmetries of the angular kernels about the equatorial axis at θ = π/2 (Dahlen & Tromp 1998):

$$\begin{array}{cc}& Y_l^{N,m\ast }={(-1)}^{-m-N}{Y}_l^{-N,-m},\hfill \end{array}$$(D.22)

$$\begin{array}{cc}& P_l^{N,-m}(\mu )={(-1)}^{l+N}{P}_l^{N,m}(-\mu ),\hfill \end{array}$$(D.23)

$$\begin{array}{cc}& P_l^{-N,m}(\mu )={(-1)}^{l+m}{P}_l^{N,m}(-\mu ).\hfill \end{array}$$(D.24)

Appendix E: Wigner 3-j symbols

The Wigner 3-j symbols provide the maximal symmetric form for expressing the addition of three angular momenta resulting in zero total momentum. They are closely related to the Clebsch-Gordon coefficients. The most symmetric form is given by (Racah 1942; Edmonds 1960):

$$\begin{array}{cc}& \left(\begin{array}{ccc}{l}_1& l_2& l_3\\ m_1& m_2& m_3\end{array}\right)\hfill \\ \hfill & ={(-1)}^{l_1-l_2-m_3}{(2l_3+1)}^{-1/2}(l_1{m}_1{l}_2{m}_2|l_1{l}_2{l}_3-m_3)\hfill \\ \hfill & ={(-1)}^{l_2+l_3+m_1}\frac{\sqrt{(l_1+l_2-l_3)!(l_1-l_2+l_3)!(-l_1+l_2+l_3)!}}{\sqrt{(l_1+l_2+l_3+1)!}}\hfill \\ \hfill & ·\sqrt{(l_1+m_1)!(l_1-m_1)!(l_2+m_2)!(l_2-m_2)!}\hfill \\ \hfill & ·\sqrt{(l_3+m_3)!(l_3-m_3)!}{\displaystyle \sum _k}\frac{{(-1)}^k}{D_k}.\hfill \end{array}$$(E.1)

Here, l_i are non-negative integers, −l_i ≤ m_i ≤ l_i are integers, and

$$\begin{array}{cc}\hfill D_k=& k!(-l_1+l_2+l_3-k)!(l_2-m_2-k)!(l_3+m_3-k)!\hfill \\ \hfill & ·(l_1-l_2-m_3+k)!(l_1-l_3+m_2+k)!,\hfill \end{array}$$(E.2)

where k ∈ ℕ₀ runs over all values for which D_k ≠ 0. The Wigner 3-j symbol remains unchanged for even permutations of columns:

$$\begin{array}{cc}\hfill \left(\begin{array}{ccc}{l}_1& l_2& l_3\\ m_1& m_2& m_3\end{array}\right)=& \left(\begin{array}{ccc}{l}_2& l_3& l_1\\ m_2& m_3& m_1\end{array}\right)\hfill \\ \hfill =& \left(\begin{array}{ccc}{l}_3& l_1& l_2\\ m_3& m_1& m_2\end{array}\right).\hfill \end{array}$$(E.3)

For odd permutations of the columns, the following holds:

$$\begin{array}{cc}\hfill \left(\begin{array}{ccc}{l}_1& l_3& l_2\\ m_1& m_3& m_2\end{array}\right)=& \left(\begin{array}{ccc}{l}_3& l_2& l_1\\ m_3& m_2& m_1\end{array}\right)\hfill \\ \hfill =& \left(\begin{array}{ccc}{l}_2& l_1& l_3\\ m_2& m_1& m_3\end{array}\right)\hfill \\ \hfill =& {(-1)}^{l_2+l_1+l_2}\left(\begin{array}{ccc}{l}_3& l_1& l_2\\ m_3& m_1& m_2\end{array}\right)\hfill \\ \hfill =& {(-1)}^{l_1+l_2+l_3}\left(\begin{array}{ccc}{l}_1& l_2& l_3\\ m_1& m_2& m_3\end{array}\right).\hfill \end{array}$$(E.4)

Changing the sign of the bottom indices yields:

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}{l}_1& l_2& l_3\\ -m_1& -m_2& -m_3\end{array}\right)={(-1)}^{l_1+l_2+l_3}\left(\begin{array}{ccc}{l}_1& l_2& l_3\\ m_1& m_2& m_3\end{array}\right).\end{array}$$(E.5)

From this property, we conclude for m₁ = m₂ = m₃ = 0 that:

$$\begin{array}{cc}& \left(\begin{array}{ccc}{l}_1& l_2& l_3\\ 0& 0& 0\end{array}\right)={(-1)}^{l_1+l_2+l_3}\left(\begin{array}{ccc}{l}_1& l_2& l_3\\ 0& 0& 0\end{array}\right)\hfill \\ \hfill \Rightarrow & \left(\begin{array}{ccc}{l}_1& l_2& l_3\\ 0& 0& 0\end{array}\right)=0\text{if}{l}_1+l_2+l_3\text{is odd.}\hfill \end{array}$$(E.6)

The coupling of angular momenta obeys certain selection rules, which are also expressed in the Wigner 3-j symbols. These rules are summarised as follows:

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}{l}_1& l_2& l_3\\ m_1& m_2& m_3\end{array}\right)\ne 0,\end{array}$$(E.7)

if and only if:

$$\begin{array}{cc}& m_1+m_2+m_3=0,\hfill \end{array}$$(E.8)

$$\begin{array}{cc}& |m_i|\le l_i\forall i,\hfill \end{array}$$(E.9)

$$\begin{array}{cc}& |l_i-l_j|\le l_k\le l_i+l_j\text{for}i,j,k=1,2,3,\hfill \end{array}$$(E.10)

$$\begin{array}{cc}& l_1+l_2+l_3\text{is}\text{(even) integer}\text{(if}{m}_1=m_2=m_3=0\text{)}.\hfill \end{array}$$(E.11)

Appendix F: Convergence of general solution for integrals over product of three GSHs

We briefly demonstrate how the general solution in equation (56) for integrals over the product of three GSHs converges to the expression presented by Dahlen & Tromp (1998), given by equation (51), in the case N₁ = N₂ + N₃.

The orthonormality relation for two generalised Legendre functions of the same order m and upper index N is given by Dahlen & Tromp (1998) as:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{-1}^1}{P}_l^{N,m}(\mu )P_{l^{\prime }}^{N,m}(\mu )\text{d}\mu =\frac{2}{2l+1}{\delta }_{l^{\prime }l}.\end{array}$$(F.1)

For N = 0, m = 0 and l = 0, this simplifies to:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{-1}^1}{P}_{l^{\prime }}^{0,0}(\mu )\text{d}\mu =2{\delta }_{l^{\prime }0},\end{array}$$(F.2)

since P₀^0, 0 = 1. This corresponds to the result of the integral over the generalised Legendre function in equation (56) for the special case N₁ = N₂ + N₃. The Kronecker delta δ_l′0 eliminates the sum over l′ in equation (56), resulting in two Wigner 3-j symbols of the form:

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}0& l_1& l\\ 0& m_1& -m_1\end{array}\right)\left(\begin{array}{ccc}0& l_1& l\\ 0& N_1& -N_1\end{array}\right).\end{array}$$

Applying the triangle condition (E.10) of the Wigner 3-j symbols restricts the second sum over l, yielding a factor δ_l1l. The two Wigner 3-j symbols then reduce to:

$$\begin{array}{cc}& \left(\begin{array}{ccc}0& l_1& l_1\\ 0& N_1& -N_1\end{array}\right)\left(\begin{array}{ccc}0& l_1& l_1\\ 0& m_1& -m_1\end{array}\right)\hfill \\ \hfill & =\left(\begin{array}{ccc}{l}_1& 0& l_1\\ -N_1& 0& N_1\end{array}\right)\left(\begin{array}{ccc}{l}_1& 0& l_1\\ -m_1& 0& m_1\end{array}\right)\hfill \\ \hfill & ={(-1)}^{N_1+m_1}\frac{(2l_1)!}{(2l_1+1)!},\hfill \end{array}$$(F.3)

where we have used equation (E.3) and equation (C.229) in Dahlen & Tromp (1998). Combining the above results recovers equation (51).

Appendix G: Angular kernels

Section 4 presents the general matrix element for a toroidal magnetic field. The terms in the general matrix element, as shown in equation (50), are sorted by their angular contributions. This yields 23 distinct angular kernels, denoted S₁ through S₂₃. These kernels include contributions from the GSHs associated with the two eigenfunctions and the magnetic field components arising from the Lorentz force. They also include trigonometric functions from spherical coordinates and derivatives of the magnetic fields and eigenfunctions in spherical geometry. Although the radial part corresponding to the angular kernel S₁₂ cancels out, we list it here for completeness. The angular kernels are defined as:

$$\begin{array}{cc}\hfill S_1& =Y_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.1)

$$\begin{array}{cc}\hfill S_2& =cot\theta Y_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \theta }{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.2)

$$\begin{array}{cc}\hfill S_3& =Y_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^2}{\partial {\theta }^2}{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.3)

$$\begin{array}{cc}\hfill S_4& =\frac{1}{{sin}^2\theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^2}{\partial {\varphi }^2}{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.4)

$$\begin{array}{cc}\hfill S_5& =Y_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \theta }{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.5)

$$\begin{array}{cc}\hfill S_6& =Y_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \theta }{Y}_l^{0,m}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.6)

$$\begin{array}{cc}\hfill S_7& =cot\theta \frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.7)

$$\begin{array}{cc}\hfill S_8& =\frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{0,m}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.8)

$$\begin{array}{cc}\hfill S_9& ={cot}^2\theta \frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \theta }{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.9)

$$\begin{array}{cc}\hfill S_{10}& =cot\theta \frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^2}{\partial {\theta }^2}{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.10)

$$\begin{array}{cc}\hfill S_{11}& =\frac{cot\theta }{{sin}^2\theta }\frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^2}{\partial {\varphi }^2}{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.11)

$$\begin{array}{cc}\hfill S_{12}& =cot\theta \frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \theta }{Y}_l^{0,m}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.12)

$$\begin{array}{cc}\hfill S_{13}& =\frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^2}{\partial {\theta }^2}{Y}_l^{0,m}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.13)

$$\begin{array}{cc}\hfill S_{14}& =\frac{1}{{sin}^2\theta }\frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^2}{\partial {\varphi }^2}{Y}_l^{0,m}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.14)

$$\begin{array}{cc}\hfill S_{15}& =cot\theta \frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \theta }{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.15)

$$\begin{array}{cc}\hfill S_{16}& =\frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \theta }{Y}_l^{0,m}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_s^{0,0}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.16)

$$\begin{array}{cc}\hfill S_{17}& =\frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \theta }{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.17)

$$\begin{array}{cc}\hfill S_{18}& =\frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^3}{\partial {\theta }^3}{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.18)

$$\begin{array}{cc}\hfill S_{19}& =\frac{1}{{sin}^2\theta }\frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^3}{\partial \theta \partial {\varphi }^2}{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.19)

$$\begin{array}{cc}\hfill S_{20}& =\frac{\partial }{\partial \theta }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \theta }{Y}_l^{0,m}\frac{\partial {}^3}{\partial {\theta }^3}{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.20)

$$\begin{array}{cc}\hfill S_{21}& =\frac{1}{{sin}^2\theta }\frac{\partial }{\partial \varphi }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial }{\partial \varphi }{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.21)

$$\begin{array}{cc}\hfill S_{22}& =\frac{cot\theta }{{sin}^2\theta }\frac{\partial }{\partial \varphi }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^2}{\partial \theta \partial \varphi }{Y}_l^{0,m}\frac{\partial }{\partial \theta }{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}\hfill \end{array}$$(G.22)

$$\begin{array}{cc}\hfill S_{23}& =\frac{1}{{sin}^2\theta }\frac{\partial }{\partial \varphi }{Y}_{l^{\prime }}^{0,m^{\prime }\ast }\frac{\partial {}^2}{\partial \theta \partial \varphi }{Y}_l^{0,m}\frac{\partial {}^2}{\partial {\theta }^2}{Y}_s^{0,0}\frac{\partial }{\partial \theta }{Y}_{s^{\prime }}^{0,0}.\hfill \end{array}$$(G.23)

In all angular kernels, we express the derivatives of the GSHs using equations (D.5) – (D.9) and eliminate the trigonometric functions using equations (D.11), (D.12), (D.14), and (D.15). This leads to angular kernels depending only on constant factors and products of GSHs. The equations derived in Section 4.1 allow us to perform the angular integral over these products of GSHs analytically. The explicit expressions of the angular kernels are listed below, incorporating the symmetry relations for the constants and the GSHs, as provided by equations (D.7) and (D.17), respectively, to improve readability:

$$\begin{array}{cc}\hfill S_1& =c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \end{array}$$(G.24)

$$\begin{array}{cc}\hfill S_2& =\frac{1}{4}{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_s^0{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}-c_s^0{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +c_s^2{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}-c_s^2{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.25)

$$\begin{array}{cc}\hfill S_3& =\frac{1}{4}{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_l^2{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{2,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}+c_l^2{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{-2,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -2c_l^0{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.26)

$$\begin{array}{cc}\hfill S_4& =-\frac{1}{4}{c}_l^0{c}_{l^{\prime }}^0{c}_s^0{c}_{s^{\prime }}^0(Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}+Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}+Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.27)

$$\begin{array}{cc}\hfill S_5& =\frac{1}{4}{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_{s^{\prime }}^0{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{0,0}-c_{s^{\prime }}^0{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_{s^{\prime }}^2{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{2,0}+c_{s^{\prime }}^2{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{2,0})\hfill \end{array}$$(G.28)

$$\begin{array}{cc}\hfill S_6& =\frac{1}{4}{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_s^0{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}-c_s^0{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -c_s^2{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}+c_s^2{Y}_{l^{\prime }}^{0,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.29)

$$\begin{array}{cc}\hfill S_7& =\frac{1}{4}{c}_{l^{\prime }}^0{c}_s^0{c}_{s^{\prime }}^0(c_s^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}+c_s^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -c_s^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}-c_s^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.30)

$$\begin{array}{cc}\hfill S_8& =\frac{1}{4}{c}_{l^{\prime }}^0{c}_s^0{c}_{s^{\prime }}^0(Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}-Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}+Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.31)

$$\begin{array}{cc}\hfill S_9& =\frac{1}{16}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}+c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}-c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}-c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}-c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}-c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}+c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}+c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0})\hfill \end{array}$$(G.32)

$$\begin{array}{cc}\hfill S_{10}& =\frac{1}{16}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_l^2{c}_s^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-2,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +c_l^2{c}_s^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-2,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}-2c_l^0{c}_s^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -2c_l^0{c}_s^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}+c_l^2{c}_s^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{2,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +c_l^2{c}_s^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{2,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}-c_l^2{c}_s^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-2,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -c_l^2{c}_s^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-2,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}+2c_l^0{c}_s^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +2c_l^0{c}_s^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}-c_l^2{c}_s^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{2,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -c_l^2{c}_s^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{2,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.33)

$$\begin{array}{cc}\hfill S_{11}& =\frac{m}{8}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(\sqrt{\frac{\pi }{3}}[c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}]\hfill \\ \hfill & -\sqrt{\frac{\pi }{6}}[c_s^0(c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0)Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2(c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0)Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0(c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0)Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2(c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0)Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0(c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0)Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2(c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0)Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0(c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0)Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2(c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0)Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}])\hfill \end{array}$$(G.34)

$$\begin{array}{cc}\hfill S_{12}& =\frac{1}{8}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}-c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}-c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}+c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}-c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}+c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}+c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}-c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0})\hfill \end{array}$$(G.35)

$$\begin{array}{cc}\hfill S_{13}& =\frac{1}{16}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_l^2{c}_s^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-2,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -c_l^2{c}_s^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-2,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}-2c_l^0{c}_s^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +2c_l^0{c}_s^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}+c_l^2{c}_s^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{2,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -c_l^2{c}_s^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{2,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}-c_l^2{c}_s^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-2,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +c_l^2{c}_s^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-2,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}+2c_l^0{c}_s^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -2c_l^0{c}_s^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{0,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}-c_l^2{c}_s^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{2,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +c_l^2{c}_s^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{2,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.36)

$$\begin{array}{cc}\hfill S_{14}& =\frac{m}{8}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(\sqrt{\frac{\pi }{3}}[c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}]\hfill \\ \hfill & -\sqrt{\frac{\pi }{6}}[c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}])\hfill \end{array}$$(G.37)

$$\begin{array}{cc}\hfill S_{15}& =\frac{1}{16}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}+c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}-c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}-c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}-c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}-c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}+c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}+c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0})\hfill \end{array}$$(G.38)

$$\begin{array}{cc}\hfill S_{16}& =\frac{1}{16}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}-c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}-c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}+c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}-c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}+c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}+c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}-c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0})\hfill \end{array}$$(G.39)

$$\begin{array}{cc}\hfill S_{17}& =\frac{1}{4}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}-Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}+Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.40)

$$\begin{array}{cc}\hfill S_{18}& =\frac{1}{16}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(c_l^2{c}_l^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-3,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -[c_l^2{c}_l^2+2c_l^0{c}_l^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +[c_l^2{c}_l^2+2c_l^0{c}_l^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -c_l^2{c}_l^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{3,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}-c_l^2{c}_l^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-3,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +[c_l^2{c}_l^2+2c_l^0{c}_l^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -[c_l^2{c}_l^2+2c_l^0{c}_l^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +c_l^2{c}_l^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{3,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.41)

$$\begin{array}{cc}\hfill S_{19}& =-\frac{m^2}{4}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(\frac{\pi }{3}[c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}]\hfill \\ \hfill & -\frac{\pi }{3\sqrt{2}}[[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_s^3{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_s^3{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_s^3{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_s^3{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_s^3{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_s^3{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_s^3{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_s^3{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_1^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}]\hfill \\ \hfill & +\frac{\pi }{6}[[c_s^2{c}_s^2-2c_s^0{c}_s^0][c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_s^3[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_s^3{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -[c_s^2{c}_s^2-2c_s^0{c}_s^0][c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_s^3[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_s^3{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -[c_s^2{c}_s^2-2c_s^0{c}_s^0][c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_s^3[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_s^3{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +[c_s^2{c}_s^2-2c_s^0{c}_s^0][c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -[c_s^2{c}_s^2-2c_s^0{c}_s^0]c_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_s^3[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_s^3{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_1^{1,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}])\hfill \end{array}$$(G.42)

$$\begin{array}{cc}\hfill S_{20}& =-\frac{1}{16}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0([c_s^2{c}_s^2+2c_s^0{c}_s^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -c_s^2{c}_s^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -[c_s^2{c}_s^2+2c_s^0{c}_s^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +c_s^2{c}_s^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -[c_s^2{c}_s^2+2c_s^0{c}_s^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +c_s^2{c}_s^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{3,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +[c_s^2{c}_s^2+2c_s^0{c}_s^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & -c_s^2{c}_s^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{3,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.43)

$$\begin{array}{cc}\hfill S_{21}& =\frac{1}{4}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}+Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}\hfill \\ \hfill & +Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0}+Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{1,0}{Y}_{s^{\prime }}^{1,0})\hfill \end{array}$$(G.44)

$$\begin{array}{cc}\hfill S_{22}& =-\frac{m}{8}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(\sqrt{\frac{\pi }{3}}[c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}]\hfill \\ \hfill & -\sqrt{\frac{\pi }{6}}[c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}])\hfill \end{array}$$(G.45)

$$\begin{array}{cc}\hfill S_{23}& =-\frac{m}{8}{c}_{l^{\prime }}^0{c}_l^0{c}_s^0{c}_{s^{\prime }}^0(\sqrt{\frac{\pi }{3}}[c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^0{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{0,0}{Y}_1^{0,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{2,0}{Y}_1^{0,0}]\hfill \\ \hfill & -\sqrt{\frac{\pi }{6}}[c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{-1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{-1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^0[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^0{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{0,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}\hfill \\ \hfill & +c_s^2[c_{s^{\prime }}^2{c}_{s^{\prime }}^2-2c_{s^{\prime }}^0{c}_{s^{\prime }}^0]Y_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{1,0}{Y}_1^{1,0}\hfill \\ \hfill & -c_s^2{c}_{s^{\prime }}^2{c}_{s^{\prime }}^3{Y}_{l^{\prime }}^{1,m^{\prime }\ast }{Y}_l^{1,m}{Y}_s^{2,0}{Y}_{s^{\prime }}^{3,0}{Y}_1^{1,0}]).\hfill \end{array}$$(G.46)

All Figures