The Direct Effect of Toroidal Magnetic Fields on Stellar Oscillations: An Analytical Expression for the General Matrix Element

Where is the solar dynamo located and what is its modus operandi? These are still open questions in solar physics. Helio- and asteroseismology can help answer them by enabling us to study solar and stellar internal structures through global oscillations. The properties of solar and stellar acoustic modes are changing with the level of magnetic activity. However, until now, the inference on subsurface magnetic fields with seismic measures has been very limited. The aim of this paper is to develop a formalism to calculate the effect of large-scale toroidal magnetic fields on solar and stellar global oscillation eigenfunctions and eigenfrequencies. If the Lorentz force is added to the equilibrium equation of motion, stellar eigenmodes can couple. In quasi-degenerate perturbation theory, this coupling, also known as the direct effect, can be quantified by the general matrix element. We present the analytical expression of the matrix element for a superposition of subsurface zonal toroidal magnetic field configurations. The matrix element is important for forward calculations of perturbed solar and stellar eigenfunctions and frequency perturbations. The results presented here will help to ascertain solar and stellar large-scale subsurface magnetic fields, and their geometric configuration, strength, and their change over the course of activity cycles.

1. INTRODUCTION So far, the main benchmarks for dynamo simulations are the correct cycle lengths as well as the properties and distribution of active regions over the simulated cycle (e.g., Charbonneau 2010Charbonneau , 2013. The strength and configuration of large-scale flow fields are important ingredients of modern simulations of solar and stellar flux transport dynamos (e.g., Choudhuri et al. 1995;Dikpati and Gilman 2009;Karak et al. 2014). Seismology has contributed to the refinement of dynamo simulations by mapping the differential rotation in the Sun (Schou et al. 1998;Howe 2009) and is beginning to do so in stars (Beck et al. 2012;Nielsen et al. 2015;Hekker and Christensen-Dalsgaard 2017). The meridional circulation is another crucial element in these dynamo models, which can be assessed by helioseismic investigations (Giles et al. 1997;Braun and Fan 1998;Haber et al. 2002;Zhao et al. 2013). If it were made possible to measure the location, geometrical configuration, and strength of the magnetic field distribution within the Sun and stars, this would add a benchmark of supreme importance for dynamo simulation.
In the Sun, a large fraction of the total magnetic flux is assumed to be stored in the subsurface toroidal component of the magnetic field (Charbonneau 2010, and references therein). The unsigned magnetic flux attributed to active regions amounts to ∼10 25 Mx over a solar cycle, while the polar cap flux totals to just ∼10 22 Mx (Charbonneau 2010). As the active regions are thought to originate from deep-seated toroidal flux concentrations (e.g., Caligari et al. 1995;Fan 2009), it is reasonable to concentrate on the toroidal field in a first step, as we will do in this article.
The effect of a magnetic field on stellar acoustic oscillations is twofold: the direct effect, or, as we shall also call it in the following, the coupling of oscillation modes, and the indirect effect. In this article, we focus on the direct effect. The indirect effect, which is due to additional forces -in comparison to the equilibrium stellar model without magnetic fields -perturbs stellar structural quantities, such as sound speed and density. Thus, the cavities where resonant acoustic modes exist are changed, which leads to changes in frequencies and eigenfunctions. The effect of a magnetic field on stellar structure was studied by, e.g., Mestel and Moss (1977), Mathis and Zahn (2005), Duez et al. (2008Duez et al. ( , 2010. We use an ansatz from quasi-degenerate perturbation theory to calculate the strength of the coupling between stellar oscillation modes (Lavely and Ritzwoller 1992), which is due to a superposition of large-scale axisymmetric zonal toroidal magnetic fields. The mode coupling leads to frequency shifts, frequency splittings, and distortions of the mode eigenfunctions. The strength of the coupling is quantified by the general matrix element, which is used to construct the supermatrix (Lavely and Ritzwoller 1992). Solving the eigenvalue problem for the supermatrix yields two essential results. First, the perturbations to the equilibrium frequencies are given by the eigenvalues. Second, the eigenvectors contain the expansion coefficients of the perturbed eigenfunction as components.
Previous studies with large-scale flows as the source of perturbation have shown that the analysis of perturbed eigenfrequencies and eigenfunctions can be employed to infer the geometrical configuration and strength of that perturbation. Ritzwoller and Lavely (1991) developed a formalism to invert global helioseismic data for differential rotation. This was expanded by Lavely and Ritzwoller (1992) to specifically include global-scale convection and expanded in such a way that general flow fields can be taken as a source of perturbation. With their theory, the perturbation of the mode eigenfrequencies and the perturbed eigenfunctions can be computed. In that framework, the effect of the meridional circulation on solar eigenmodes was studied by, e.g., Roth and Stix (2008). Later, Schad et al. (2011) advanced the method by introducing the amplitude ratio of oscillations in the Fourier domain as a helioseismic measurement quantity. With this,  inferred a multicellular meridional circulation in both depth and latitude by analyzing perturbed solar eigenfunctions.
Using a different perturbation approach, Gough and Thompson (1990) presented a formalism to calculate the eigenfrequencies of a stellar model with a magnetic field and rotation as perturbations. This formalism included both the direct and indirect effects of the magnetic field. Antia et al. (2000) refined their approach by also including effects from general relativity and second-order effects of rotation. Their main results were forward calculations of helioseismic splitting coefficients. They could limit the strength of a toroidal magnetic field near the Sun's tachocline to 300 kG. Baldner et al. (2009) used the theory developed by Gough and Thompson (1990) and Antia et al. (2000) to probe the solar interior for magnetic field concentrations by comparing the calculated splitting coefficients with the observed values. Their results included two very shallow toroidal field distributions and a poloidal field. However, due to the shallowness of the inferred fields, they could not draw conclusions regarding dynamo models, which mostly rely on deeper seated field concentrations (e.g., Charbonneau 2014). Recently, Hanasoge (2017) developed a formalism to calculate the sensitivity functions of seismic measurements to the Lorentz stress tensor in the stellar interior.
A seismic investigation of the effect of magnetic fields on the oscillations in the Sun and solar-like stars within the perturbation approach developed by Lavely and Ritzwoller (1992) has not been done so far. The first step of extending their analytical groundwork to fully include magnetic fields is now presented in this paper. We begin by deriving the perturbed equation of motion for a non-rotating, non-magnetic star in Section 2. The necessary essentials of quasidegenerate perturbation theory are reviewed in Section 3. Section 4 is dedicated to the derivation of the general matrix element for a magnetic field, before we give its explicit analytical expression of it in spherical coordinates in Section 5. We discuss our findings and conclude in Section 6.

THE EQUILIBRIUM MODEL
In this Section, we summarize the derivation of the perturbed equation of motion, because it will help us to trace our steps in later Sections. A detailed treatment of this can be found in, e.g., Unno et al. (1989) and Aerts et al. (2010).
The equation of motion of a parcel of gas in a non-rotating, static, and spherically symmetric star without a magnetic field reads where v is the velocity due to displacements from the equilibrium position, p is the pressure, g is the gravitational acceleration, and ρ is the mass density. For such a star, a standard model can be computed, e.g. for the Sun Model S by Christensen-Dalsgaard et al. (1996). The model supplies the eigenfunctions and eigenfrequencies of adiabatic oscillations as a solution of the eigenvalue equation, which is briefly derived here. The Lagrangian perturbation ∆ of a quantity Q is given by where δQ is the Eulerian variation of Q and (ξ · ∇)Q is the term due to the advection of Q by the displacement ξ.
The Lagrangian perturbation commutes with the total time derivative: Furthermore, the following identifications of Lagrangian perturbations hold, which we need to rewrite the perturbed equation of motion: where r is the position vector. Here, the material time derivative is equal to the partial time derivative since we consider a star without any velocity fields, such as rotation, convection, or meridional circulation. We start by dividing Equation (1) by ρ, taking the Lagrangian perturbation ∆ of the whole equation, and multiplying by the unperturbed density ρ 0 . From now on, unperturbed quantities are labeled with a sub-or superscript 0. These operations, together with Equations (3)-(6), yield A change between the gravitational potential and gravitational acceleration can be done via −∇Φ = g. By using a modal ansatz for the displacements, and making use of the identification ρ 0 g 0 = ∇p 0 , the perturbed equation of motion finally reads We identify the right-hand side with the operator L 0 (ξ 0 ): This equation describes the response of the stellar model to a displacement ξ of a parcel of gas given by the Lagrangian perturbation of Equation (1). This is the eigenvalue equation we sought to derive. The superscript 0 on the eigenfunctions is added here to signify that they are eigenfunctions of the equilibrium model. By adding additional forces or velocity fields to the equation of motion, we can probe the response of the model and its eigenfunctions to this perturbation. It may be noted that Equation (10) is equivalent to Equation B15 in Lavely and Ritzwoller (1992) and to Equation 28 in Lynden-Bell and Ostriker (1967).
In this paper, all calculations are performed in spherical polar coordinates. In the following, r stands for the radial distance from the origin, θ signifies the colatitude, and φ is the azimuthal angle. As we are working with a spherically symmetric star, the r, θ, and φ components of the eigenfunctions can be written as Here, ξ r (r) and ξ h (r) are the real-valued scalar radial and horizontal displacement amplitudes, respectively. The generalized spherical harmonic functions Y N,m l are described in detail in Appendix D. The set of eigenmodes {ξ k }, where k = (l, m, n), with eigenfrequencies ω k , gives the solutions to Equation (11). The integer l is the harmonic degree, m is the azimuthal order, with |m| ≤ l, and n is the radial order of the eigenfunction. As we are considering a non-rotating star, modes of the same radial order and harmonic degree (n, l) but of different azimuthal order m have the same frequency and form a degenerate multiplet (see, e.g., Unno et al. 1989): ω k = ω nl .
The unperturbed eigenfunctions are orthogonal, and normalized with the normalization constant The integral in Equation (13) extends over the whole stellar volume and the integral in Equation (14) extends over the stellar radius.
3. ESSENTIALS OF QUASI-DEGENERATE PERTURBATION THEORY We briefly review the required essentials of quasi-degenerate perturbation theory here. More detailed treatments can be found in, e.g., Lavely and Ritzwoller (1992), Roth (2002), or . We introduce the following perturbation expansions into Equation (11): They express the effect on the eigenfrequencies, eigenfunctions, and equilibrium equation of motion by the operator L 1 that accounts for departures from the equilibrium model. The auxiliary quantity signifies a small perturbation and helps to keep track of the order of perturbation. The squared eigenfrequency perturbation in first order is given by (ω 2 1 ) j . The perturbation of the equation of motion leads to a coupling of eigenmodes. The coupling due to, e.g., poloidal flow cells was studied by Herzberg (2016) and that due to a meridional flow in the Sun by Schad et al. (2011). In the framework of quasi-degenerate perturbation theory, the perturbed eigenfunctionsξ j to zeroth and first orders, i.e.,ξ 0 j andξ 1 j , can be described as a weighted sum over the unperturbed eigenfunctions ξ 0 k , for which the modes k are within a coupling set K or K ⊥ , multiplied by a coupling coefficient: Which modes are within the coupling set K or its complement K ⊥ depends on the geometrical configuration of the perturbation, which entails angular momentum selection rules, on the considered modes, as well as on the allowed difference in the frequency of the modes. When, e.g., an axisymmetric meridional poloidal flow is the source of perturbation, only modes of the same azimuthal order can couple due to selection rules (e.g., Lavely and Ritzwoller 1992;Roth 2002). We will discuss the selection rules for toroidal magnetic fields in Section 5.1. The coupling coefficients c jk and b jk may be complex valued. The frequency difference of modes within K is restricted by the quasi-degeneracy condition: Here, ∆ω 2 sets the frequency range of the coupling modes within K. The reference frequency ω ref is typically chosen to be equal or close to the central frequency of a multiplet j. Increasing ∆ω 2 results in more modes within the coupling set K and therefore in a more accurate expression for the perturbed quantity. This has to be traded off against the higher computational effort necessary for the computation. We can focus on the determination of the coefficients c jk , which are needed to calculate the perturbed eigenfunction in zeroth orderξ 0 j in Equation (16a). With the perturbation expansion, Equations (15a)-(15c), the following eigenvalue equation can be found (Lavely and Ritzwoller 1992): Here, it was used that a mode k is either in K or K ⊥ , i.e., the equations for the coefficients c jk and b jk decouple. Now, the sought-for coefficients c jk appear as the elements of the eigenvectors of the matrix Z. This matrix is called supermatrix, and its elements are determined from It can be shown that the general matrix element H k k is given by where L 1 is the operator of the perturbing force that causes the coupling, introduced in Equation (15a). The eigenvectors and eigenvalues of the supermatrix completely determine the perturbed eigenfunctions and frequencies. The squared frequency correction of a mode k is given by the eigenvalue, which belongs exactly to the eigenvector of Z that holds the expansion coefficients c jk k∈K for the construction of the perturbed eigenfunction of the mode.

THE GENERAL MATRIX ELEMENT FOR A MAGNETIC FIELD
As a perturbation to the equilibrium, we add the Lorentz force F L = 1 c j × B, where the current is given by to the right-hand side of the equation of motion, Equation (1). We use the CGS-Gauss system for all calculations. Following the same steps as in Section 2 leads to where we made the identification This expression is the desired perturbation operator for a magnetic field, which can now be inserted into Equation (20) with L B (ξ 0 ) = L 1 (ξ 0 ). Doing so yields the general matrix element for the Lorentz force exerted by a magnetic field of general configuration: where Equation (C1) was used. The Eulerian perturbation of the magnetic field δB k can be expressed as The subscript k of the eigenmodes signifies that they are solutions to Equation (11). This is in contrast to ξ in the derivation of identity (25) -see Appendix A -which can be any displacement. Identity (25) is equivalent to the induction equation of ideal magnetohydrodynamics for small perturbations in the magnetic field and the velocity field. The unperturbed solar or stellar model has a (2l + 1)-fold degeneracy in each frequency multiplet (n, l). This degeneracy can easily be lifted by introducing a differential rotation profile as a perturbation to the equilibrium model, and then applying the theory of Lavely and Ritzwoller (1992) as discussed in, e.g., Roth (2002).

Modeling the Toroidal Magnetic Field
So far, B has been a magnetic field of any configuration. We now restrict B to a superposition of purely toroidal fields without a φ dependence, which is represented in generalized spherical harmonics as Here, a s (r) are the radial profiles of the toroidal field components and Y 0,0 s (θ, φ) are the generalized spherical harmonic functions. Some useful properties of these functions are collected in Appendix D. The index range of s extends over those even positive integers, which are contributing to the magnetic field of the model. The restriction to even integers ensures opposite polarities in the northern and southern hemispheres.

The Matrix Element of a Toroidal Magnetic Field
With the specialization to a superposition of toroidal magnetic field components (Equation (26)), the general matrix element from Equation (24) is now given by where s and s extend over the same index range. Due to the distributivity of the vector product and the dot product as well as due to the sum rule for the integral, this simplifies to where we applied Equation (25). We will now calculate the extended form of this equation in spherical polar coordinates.
5. THE ANALYTICAL FORM OF THE GENERAL MATRIX ELEMENT To simplify notation and avoid indices, we consider one summand from Equation (28) for fixed s, s . We omit the indices s, s and use where a s = a and a s =â. We use the same notation for the eigenmodes, which appear twice in Equation (28), once due to the Lagrangian perturbation of the Lorentz force and once for the eigenmode that is multiplied from the left in Equation (28). The eigenmode that is multiplied from the left is assigned the hat, i.e., its radial component is written . We line out the necessary calculations in Appendix C. As a result, we find the following expression for the general matrix element given by Equation (28): The factors with a radial dependence have been separated from those factors with an angular dependence. The angular parts are contained in the symbols T i with i = 1, . . . , 25. They are listed in Appendix G. The integral in Equation (31) extends over the whole stellar volume. Integration over the solid angle can be carried out by utilizing Wigner 3j symbols, as demonstrated in Appendix (F). The properties of the Wigner 3j symbols (see Appendix (E)), can be used to reduce the number of terms in the angular kernels (G1)-(G25), which leads to Equations (G27)-(G51). The summation over the components of a superposition of toroidal magnetic fields as in Equation (25) can now easily be applied again by assigning the indices s, s to the two magnetic fields as defined in Equations (29)-(30).

Selection Rules
From the properties of the integral over four generalized spherical harmonics, given by Equation (F3), we find the following selection rules for mode coupling: two modes can only couple if their harmonic degrees satisfy This follows from property (E6c) of the Wigner 3j symbol. The azimuthal order of the coupling modes must satisfy which is a consequence of the property (E6a) of the Wigner 3j symbol and the fact that the azimuthal order was set to zero for the toroidal magnetic field.
6. DISCUSSION AND CONCLUSION Magnetic fields in stars lead to the coupling of stellar oscillation eigenmodes. This coupling perturbs the mode frequencies and eigenfunctions of these acoustic oscillations. With an ansatz from quasi-degenerate perturbation theory, in which the magnetic field is treated as a small perturbation to the equilibrium stellar model, we presented a detailed derivation of the general matrix element. In this work, we focused on the direct effect caused by a superposition of toroidal magnetic fields. The direct effect describes the mode coupling due to the Lorentz force, whereas the indirect effect, which we did not consider in this article, is due to the perturbation of the mode cavities by the Lorentz force.
Both the direct and indirect effects have to be accounted for, if a full description of the perturbed mode frequencies and eigenfunctions is to be reached. The combination of both effects will be considered in an upcoming study.
The analytical expression of the general matrix element was obtained with the use of generalized spherical harmonics and the calculation of the angular integral over their product. The resulting general matrix element, presented in Equation (31), can be used to carry out forward calculations of the direct effect of arbitrary subsurface zonal toroidal field configurations on the stellar eigenmodes. From the general matrix elements, the supermatrix can be constructed. Its eigenvalues and eigenvectors hold the information on the mode frequency perturbations and the perturbed eigenfunctions, respectively.
From an inspection of the general matrix element and the selection rules for coupling of angular momenta, we find that toroidal magnetic fields couple only modes that are of equal azimuthal order (Equation (33)). Also, only modes whose harmonic degrees maximally differ by the sum of the harmonic degrees of the toroidal magnetic field configuration can couple (Equation (32)). It is noteworthy that the general matrix element and hence the supermatrix Z is not Hermitian. This is in stark contrast to the supermatrix for toroidal flows, e.g., rotation (Lavely and Ritzwoller 1992), or poloidal flows as the meridional circulation (Roth and Stix 2008), where Z is Hermitian. The Hermiticity of the supermatrix for flows ensures that its eigenvalues and consequently the frequency perturbations are purely real. This need not be the case if the perturbation is a toroidal magnetic field. The non-Hermiticity of the supermatrix, or more precisely its asymmetry as it is completely real-valued in the case of magnetic fields, allows complex-valued eigenvalues and thus complex frequency perturbations. The asymmetry can be easily seen in Equation (31), considering, for example, the radially dependent part of the first term, which includes the angular factor T 1 . Clearly, the terms are not symmetrical under the exchange of the modes k and k, which renders the general matrix element asymmetrical. A complex-valued frequency perturbation is tantamount to a frequency shift and an additional damping factor for the affected mode.
The presented matrix element can be used to study the direct effect of toroidal field configurations motivated by the simulation of solar and stellar dynamos, e.g., by Miesch and Teweldebirhan (2016), in forward calculations for their effect on the solar and stellar eigenmodes. The resulting frequency and damping perturbations can then be compared with observational data for frequency shifts (e.g., Woodard and Noyes 1985;Libbrecht and Woodard 1990;Jimenez-Reyes et al. 1998;Broomhall 2017) and mode damping rates (e.g., Jefferies et al. 1990;Komm et al. 2000;Broomhall et al. 2015) over the solar cycle as well as the frequency shifts caused by changes in stellar magnetic activity (García et al. 2010;Salabert et al. 2016;Kiefer et al. 2017). This will provide a novel tool to calibrate and test dynamo simulations and may even help magnetic field concentrations to be locate within the solar and stellar interiors. With the formalism presented here, it will also be possible to produce artificial splitting coefficients for toroidal fields. These can then also be compared to observed splitting coefficients, which may additionally help in locating subsurface magnetic field structures.
It has been shown in previous studies that the analysis of perturbed eigenfunctions can be used to infer the profile of the solar meridional flow . The results from the present work are a first step toward such an inversion with the aim to infer the solar magnetic field from perturbed eigenfunctions analogous to the work by . In a next step, also the indirect effect of the magnetic field also has to be accounted for, which is a work in progress.
In contrast to earlier studies, which aimed to determine the effect of subsurface magnetic fields in the Sun on global oscillations (e.g., Gough and Thompson 1990;Antia et al. 2000;Baldner et al. 2009), we present an analytical expression for the direct effect. In this way, it can readily be investigated how different modes are affected by the input magnetic field depending on the amplitudes of their radial and horizontal displacements, ξ r and ξ h ; their harmonic degree; frequency; and other parameters. Modes of different frequencies and harmonic degrees probe different depths of the star (e.g., Aerts et al. 2010). Depending on this lower turning point, modes are variably susceptible to perturbations by a magnetic field of a certain location and configuration. In an upcoming study, we will exploit these properties of the presented formalism by inputting different magnetic field configurations and learn how different modes are affected by different magnetic fields.
We note that the analytical expression of the general matrix element for the direct effect of a superposition of poloidal fields or general magnetic field configuration will lead to even more extensive matrix elements than the one presented here. Nonetheless, this will be a worthwhile effort and should be carried out in the future.
With the result of this paper, we contribute to a novel approach to measure subsurface magnetic field concentrations in the Sun and stars. The systematic forward modeling of frequency shifts, frequency splittings, and the distortion of mode eigenfunctions will help us to learn about the location of solar and stellar dynamos and the determination of their mode of operation.
The authors wish to thank the anonymous referee and our institute's internal referee Oskar Steiner for taking the time to review this article.
The research leading to these results received funding from the European Research Council under the European where B is the magnetic field and u is the velocity field. We introduce the perturbations where δB and v are the Eulerian perturbation to the magnetic field and the velocity field, respectively. The unperturbed magnetic field B 0 (r) is assumed to be static. We neglect large-scale velocity fields here, as we are focusing on the influence of the magnetic field. Hence, the background velocity field u 0 (r) is set to zero. The Eulerian perturbation of the velocity field v(r, t) is described as a temporal change in the displacement ξ: v = ∂ξ ∂t . (A4) Applying Equations (A2)-(A4) to Equation (A1), retaining only terms that are linear in the perturbations, and integrating over time yield the Eulerian perturbation of the magnetic field δB: B. MATHEMATICAL SUPPLEMENTS In this appendix, we list the vector identities and formulae that were used in the calculation of the general matrix element in Section 5 and Appendix C, and help the traceability of our derivations. All calculations are done in spherical polar coordinates. In the following, r denotes the radial distance from the origin, θ signifies the colatitude, and φ is the azimuthal angle.
Let A (r, θ, φ) and B (r, θ, φ) denote general vector fields. A useful identity for the curl of the cross-product of two vector fields A and B is given by The divergence of a vector field A is and its curl is The material derivative of the vector fields A and B is given by C. CALCULATION OF THE GENERAL MATRIX ELEMENT In this appendix, we derive the explicit expression for the general matrix element in spherical coordinates. We start from the Lagrangian perturbation of the Lorentz force: where Equations (2) and (21), in which the Lagrangian perturbation ∆ and the electric current j are defined, have been used, and Equation (A5) for the Eulerian perturbation of the magnetic field δB holds. Before applying the Lagrangian perturbation to the Lorentz force, we divided by ρ only to multiply by ρ 0 afterwards. Due to the Lagrangian perturbation of 1/ρ, a term with a factor of −∆ρ/ρ 0 appears. This term was rewritten by virtue of the Lagrangian perturbation of the density ∆ρ = −ρ 0 ∇ · ξ, which can be derived by linearizing the perturbed equation of mass conservation. Furthermore, c is the speed of light and ξ is an eigenfunction. We restrict the magnetic field to be of toroidal shape but do not yet specify the exact configuration: By applying Equations (A5), (B1)-(B4), and (C2) to Equation (C1), we arrive at In this step, we kept the four terms from Equation (C1) separated and marked them off with curly brackets. We also applied the notation for the magnetic field from Section 5. The magnetic field components B φ andB φ are replaced by their representation and the components of the eigenfunctions are expressed as given by Equation (12).
In the next step, the eigenfunctionξ l ,m (r, θ, φ), where the bar denotes complex conjugation, is multiplied from the left, and the integral over the whole stellar volume is applied. As the radial and horizontal displacement amplitudes are purely real, the complex conjugation on them can be dropped:ξ r (r) =ξ r (r),ξ h (r) =ξ h (r). Then, terms that carry only a radial dependence are separated from those with an angular dependence. Terms with identical angular dependence can be factored out and summarized. After all of these steps, we finally arrive at Equation (31). The factors with an angular dependence, T i , with i = 1, . . . , 25, which appear in that equation, are listed in Appendix G.

D. THE GENERALIZED SPHERICAL HARMONICS AND THEIR PROPERTIES
The generalized spherical harmonics can be used to expand tensor fields of any order. Here, they are employed in the representation of the eigenfunctions and the toroidal magnetic fields as given in Equations (12), (29), and (30). We follow the conventions for the generalized spherical harmonics of Dahlen and Tromp (1998), who largely build upon the notation and sign conventions from Phinney and Burridge (1973).
The scalar generalized spherical harmonics are defined by where γ l = (2l + 1) /4π and P N lm (cos θ) are the generalized Legendre functions. In this convention, Y m l = Y 0,m l holds, where Y m l is an ordinary spherical harmonic function of degree l and azimuthal order m. We collect all identities and relations of the generalized spherical harmonics, which are needed in the calculation of the matrix element H k k , especially in the calculation of volume integrals over products of these functions; see Appendix F and G.
The generalized spherical harmonics satisfy the recursion relation with The colatitudinal derivative of Y N,m l can be expressed in terms of Y N ±1,m From Equation (D1), it can be seen that the azimuthal derivative is given by Higher colatitudinal and azimuthal derivatives can easily be constructed by an iterative application of Equation (D4)-(D5). For example, the second colatitudinal derivative of Y N,m l is given by The Two more useful relations can be obtained by taking the colatitudinal derivative of the recursion relation (D2), which yields By setting N = 0 in Equation (D8) and using Y 1,0 1 = γ 1 P 1,0 By setting m = 0 in Equation (D8) and using Y 0,0 E. THE WIGNER 3J SYMBOL The Wigner 3j symbols are maximally symmetrized representations of the Clebsch-Gordon coefficients. They represent the coupling of three angular momenta to a resultant zero momentum (Racah 1942;Edmonds 1960). The most symmetrical form to calculate them is l 1 l 2 l 3 m 1 m 2 m 3 = (−1) l 2 +l 3 +m 1 (l 1 + l 2 − l 3 )! (l 1 − l 2 + l 3 )! (−l 1 + l 2 + l 3 )! (l 1 + l 2 + l 3 + 1)!
where k ∈ N 0 runs over all values for which D k = 0 with In general, l i and m i can be integers or half-integers. In this work, l i are always non-negative integers and m i are integers in the range −l i ≤ m i ≤ l i . We list only those properties of the Wigner 3j symbols that are useful in the calculation of the general matrix element H k k . More extensive lists can be found in, e.g., Dahlen and Tromp (1998), Edmonds (1960), and Regge (1958).
The numerical value of a 3j symbol is unchanged under an even permutation of its columns: An odd permutation of its columns is equivalent to multiplication by (−1) l 1 +l 2 +l 3 : Selection rules for the allowed ranges of the l i can be constructed from the following properties: and |m i | ≤ l i for i = 1, 2, 3, and l i − l j ≤ l k for i, j, k = 1, 2, 3.
F. INTEGRAL OF THE PRODUCTS OF GENERALIZED SPHERICAL HARMONICS ON THE UNIT SPHERE To calculate the analytical expression of the general matrix element, it is necessary to integrate the product of four, five, and six generalized spherical harmonics. In Dahlen and Tromp (1998), it is shown that the integral of a product of three generalized spherical harmonics on the unit sphere can be written as where γ k = (2k + 1) /4π. The product of two generalized spherical harmonics is given by We can use Equations (F1)-(F2) to calculate the integral over four generalized spherical harmonics: The sum extends over all values of l, N , and m, which are allowed according to Equations (E6a)-(E6c). By employing Equation (E6a), we can deduce two useful conditions for the indices of the four generalized spherical harmonics in Equation (F3) if the integral shall not vanish: By iterative application of Equations (F1)-(F3), also integrals over products of five, six, or any number of generalized spherical harmonics can also be calculated. Conditions (F4) and (F5) also hold for integrals over products with more than four generalized spherical harmonics.

G. THE ANGULAR KERNELS FOR TOROIDAL MAGNETIC FIELDS
In the general matrix element, given by Equation (31), we factored out the parts that hold an angular dependence. These factors are given by the generalized spherical harmonic functions from the two eigenfunctions and the two magnetic field components, which appear in the Lorentz force. They also contain trigonometric functions, which appear due to the derivatives in spherical geometry, and derivatives acting on the generalized spherical harmonics. In this list, we include the factors T 13 -T 16 and T 21 . In our derivations, they appear for individual terms during the steps outlined in Appendix C but cancel along the way to the final result. We chose to list them here, because they may be of interest for studies that focus on the Eulerian perturbation of the Lorentz force or if poloidal magnetic fields are included. The individual factors are given by As the integral in Equation (31) extends over the whole volume and we were able to separate the radial from the angular parts, we can carry out the radial and angular integrals individually. In the following, we give the explicit expressions for the angular kernels (G1)-(G25). In the derivation of Equations (G27)-(G51), Equations (D2)-(D6) and (D10) are used to eliminate the derivatives of the generalized spherical harmonics and the trigonometric functions. We utilize selection rule (F4) to eliminate terms whose angular integral is equal to zero for any combination of harmonic degrees and azimuthal orders of the eigenfunctions and the magnetic field components. To ease readability, we use that which can be seen from Equation (D3). Also, we set m = m, which is required by Equation (F5). With all these, we find after considerable but straightforward algebra