The Effect of Toroidal Magnetic Fields on Solar Oscillation Frequencies

Solar oscillation frequencies change with the level of magnetic activity. Localizing subsurface magnetic field concentrations in the Sun with helioseismology will help us to understand the solar dynamo. Because the magnetic fields are not considered in standard solar models, adding them to the basic equations of stellar structure changes the eigenfunctions and eigenfrequencies. We use quasi-degenerate perturbation theory to calculate the effect of toroidal magnetic fields on solar oscillation mean multiplet frequencies for six field configurations. In our calculations, we consider both the direct effect of the magnetic field, which describes the coupling of modes, and the indirect effect, which accounts for changes in stellar structure due to the magnetic field. We limit our calculations to self-coupling of modes. We find that the magnetic field affects the multiplet frequencies in a way that depends on the location and the geometry of the field inside the Sun. Comparing our theoretical results with observed shifts, we find that strong tachocline fields cannot be responsible for the observed frequency shifts of p modes over the solar cycle. We also find that part of the surface effect in helioseismic oscillation frequencies might be attributed to magnetic fields in the outer layers of the Sun. The theory presented here is also applicable to models of solar-like stars and their oscillation frequencies.


INTRODUCTION
The mapping of subsurface magnetic field concentrations and with it the pinpointing of the region in the Sun, where its dynamo operates can surely be regarded as one of the outstanding open issues of solar physics. In this article, we develop the theory for the forward calculation of the effect of a superposition of zonal toroidal magnetic fields on solar oscillation frequencies in the framework of quasi-degenerate perturbation theory. The toroidal component of the global solar magnetic field is assumed to be responsible for the bulk of phenomena associated with magnetic activity (see Charbonneau 2010;Fan 2009;Hathaway 2015, and references therein). In simulations of flux-transport solar dynamos, the energy that is stored in the toroidal component of the large-scale magnetic field is orders of magnitude larger than the energy in the poloidal field (e.g., Miesch and Teweldebirhan 2016). Hence, the restriction to purely toroidal magnetic field configurations in the present work is adequate.
From observations, it is well known that solar p-mode frequencies vary in phase with the solar activity cycle. Woodard and Noyes (1985) measured these changes in frequencies for oscillation of low harmonic degree. Later, the frequency shifts over the solar cycle were confirmed and thoroughly investigated by, e.g., Libbrecht and Woodard (1990), Jiménez-Reyes et al. (1998), and Broomhall (2017). In the Sun, these shifts are larger for modes of higher frequency (Jiménez-Reyes et al. 1998). As, e.g., Basu et al. (2012) showed, modes of higher frequency have their maximum sensitivity to structural changes in the solar interior in shallower layers than modes of low frequency. This can be used to study the change of subsurface solar activity as a function of time and depth (Basu et al. 2012).
It has been shown that the frequencies of acoustic oscillations of solar-like stars undergo changes similar to those observed on the Sun. The first time that such changes were reported for a star other than the Sun was done by García et al. (2010). They also showed that these changes are correlated with stellar magnetic activity. Since then this behavior has been observed for several more stars by Salabert et al. (2016) and Kiefer et al. (2017a). Gaining insight about the underlying magnetic field causing these changes in stars would supplement the simulation of solar and stellar dynamos greatly.
Until now, the attempts to determine the magnetic field in the Sun with helioseismology have been rather sparse. Gough and Thompson (1990) calculated multiplet shifts for low degree modes for an axisymmetric buried magnetic field in a perturbational approach. Building on their framework, Antia et al. (2000) analyzed splitting coefficients to probe the solar acoustic asphericity and magnetic fields in the convection zone. They were able to limit the magnetic field strength at the base of the convection zone to 300 kG. Baldner et al. (2009) extended this work and matched simulated splitting coefficients to their observed counterparts. They found that a superposition of two very shallow toroidal magnetic fields in the upper 1% of the convection zone and a poloidal component could explain the observed even splitting coefficients best. Dziembowski and Goode (2005) analyzed the behavior of centroid multiplet frequencies between the minimum and maximum of cycle 23. They found that the frequency increase can be explained by a less than 2% decrease in the radial component of the turbulent velocity in the convection zone. They explain the p-mode frequency increase with rising levels of activity by the inhibiting effect of the magnetic field on convection. Recently, Hanasoge (2017) developed a formalism to calculate sensitivity kernels of mode coupling to Lorentz stresses in the Sun.
In this article, we describe our forward calculations of the effect of subsurface toroidal magnetic fields on solar oscillation multiplet frequencies. For this, we use an ansatz from quasi-degenerate perturbation theory to calculate the strength of the coupling between solar oscillation modes (Lavely and Ritzwoller 1992). The coupling of initially independent modes leads to changes in mode frequencies and eigenfunctions. The distortions of the solar mode eigenfunctions by flows have previously been exploited by Schad et al. (2013) to infer a double cell profile of the meridional circulation. In the ansatz, we use here, the total perturbation by the magnetic field can be separated in a direct and an indirect contribution. The direct contribution couples the modes and thus changes their frequencies and eigenfunctions. An analytical derivation of the general matrix element of this direct contribution caused by a superposition of zonal toroidal magnetic field was recently presented by Kiefer et al. (2017b). The indirect contribution is due to the perturbation to the equilibrium stellar structural quantities, which is caused by the magnetic field. The effect of a magnetic field on stellar structure was studied by, e.g., Duez et al. (2008Duez et al. ( , 2010; Mathis and Zahn (2005); Mestel and Moss (1977).
We start by introducing the necessary theoretical background of quasi-degenerate perturbation theory in Section 2. The general matrix element for the indirect effect of toroidal magnetic fields is then derived in Section 3. We present the six magnetic field configurations we tested and the resulting multiplet shifts in Section 4. These results are discussed in Section 5 and we end the paper with our conclusions in Section 6. We include mathematical supplements in Appendix A, a detailed derivation of the projection of the Lorentz force onto spherical harmonics in Appendix B, the sensitivity kernels in Appendix C, derivations of the perturbations in stellar structural quantities in Appendix D, and additional figures for the modeled magnetic fields and the resulting frequency shifts in Appendix E.

PERTURBATION THEORY
The equations that describe the equilibrium state of a star -without flows, rotation, or magnetic field -can be solved to give a system of eigenfunctions (Christensen-Dalsgaard 2008). These eigenfunctions, with their respective eigenfrequencies, are perturbed when flows (Lavely and Ritzwoller 1992), rotation, or a magnetic field (e.g., Gough and Thompson 1990) are added to the star. If the perturbation is small enough, the perturbed eigenfunctions and eigenfrequencies can be obtained from the unperturbed eigenfunctions and eigenfrequencies with techniques from standard perturbation theory. As the spectrum of the eigenfrequencies of a solar-like star is dense, the techniques from quasi-degenerate perturbation theory can be adopted. The solutions to the eigenvalue problem define the perturbed eigenvectors and eigenfrequencies of the system. Here, Z is called the supermatrix with entries where H k k is the general matrix element, which is discussed in detail below, and the indices k = (n, l, m), k = n , l , m define the considered eigenmodes with radial orders n, n , harmonic degrees l, l , and azimuthal orders m, m . The coupling set K is made up of those modes, which satisfy the following two conditions: Firstly, their frequencies obey the quasi-degeneracy condition where the reference frequency ω ref is typically chosen equal or close to the central frequency of a multiplet k, ω k is the frequency of the considered mode, and ∆ω 2 is the width of this range. Second, the geometry of the modes, determined by their harmonic degree and azimuthal order, complies with angular momentum selection rules imposed by the configuration of the perturbation. These selection rules are given in Section 5.1 of Kiefer et al. (2017b) for the magnetic field configurations considered in this work. The eigenvector C holds the expansion coefficients of the perturbed eigenfunction ξ j : where c jk is the k'th component of C and ξ 0 k is an unperturbed eigenfunction. The matrix Λ is a diagonal matrix with the frequency perturbations δω 2 k as entries. Detailed discussions of quasi-degenerate perturbation theory with application to solar and stellar problems can be found in, e.g., Lavely and Ritzwoller (1992), Roth (2002), Schad (2013), Herzberg (2016), and Herzberg and Roth (2018).
The complete general matrix element H k k is calculated as where D k k is the general matrix element, which accounts for the mode coupling due to the magnetic field. An analytical expression of D k k for a superposition of zonal toroidal magnetic fields was recently derived by Kiefer et al. (2017b). It is given by where B s , B s are toroidal magnetic fields with harmonic degrees s and s , ξ k and ξ k are eigenfunctions, and the integral extends over the stellar volume V . The sum in Equation (6) extends over the values of harmonic degrees s and s of the investigated magnetic field model, see Equation (12). The cross-terms between configurations of unequal harmonic degree are thus included in the calculation of the matrix element. The full analytical expression of D k k is given in Equations (31), and (G62)-(G86) of Kiefer et al. (2017b). I k k in Equation (5) accounts for the modal interactions due to the perturbations in stellar structural quantities, which arise as a result of the magnetic field. We derive I k k for the case of a superposition of zonal toroidal magnetic fields in Section 3. The shift in angular frequency of a mode can be approximated by where δω 2 k is the perturbation of the squared angular frequency: which is expanded up to second order (Schad et al. 2011). In this article, we do not consider perturbations to the eigenfunctions but focus on the more accessible quantity -the shifts in mode frequency. However, the perturbation of the eigenfunctions by a superposition of zonal toroidal magnetic fields can be calculated with the theory presented here: To second order, the perturbed eigenfunction ξ k is given by where ξ 0 k is the unperturbed eigenfunction. Equations (8) and (9) are basic results of non-or quasidegenerate perturbation theory, see, e.g., Sakurai and Napolitano (2014).
In the present work, we concentrate on the self-coupling of modes. This reduces the computational effort to obtain H k k by a factor of |K| 2 , i.e., from the square of the number of coupling modes to only one general matrix element H kk . This approximation is justified, as the coupling strength decreases with increasing frequency difference and the radial and horizontal eigenfunction become less similar for modes of different radial order and harmonic degree. We expect the error introduced by this approximation to be of the order of a few percent. Schad (2013) found that for differential rotation, the self-coupling matrix elements are typically two orders of magnitude larger than the crosscoupling matrix elements. We expect this to be similar for a magnetic field as the perturbing agent. The eigenfunctions are not perturbed in the self-coupling limit, as can be seen from Equation (9). In this approximation, we find from Equations (7) and (8) that the shift in mode frequency can be calculated by

THE INDIRECT EFFECT
If a magnetic field is present, the Lorentz force has to be added to the equation of motion. Hence, the structure of the star is slightly changed compared to the equilibrium. We treat this change as a small perturbation to the nonmagnetic star. The Lorentz force for a superposition of axisymmetric zonal toroidal magnetic fields is given by with the toroidal magnetic field The second magnetic fieldB tor is introduced to keep track of the individual contributions in a superposition of components of distinct harmonic degree. Here, a(r) andâ(r) are the radial profiles of the component of the magnetic field with harmonic degrees s and s , respectively. Y 0 s (θ) is a spherical harmonic function of degree s and azimuthal order 0. In this article, we will be referring to relations and properties of the spherical harmonic functions and the generalized spherical harmonic functions Y N,m l , which can be found in Appendix (D) of Kiefer et al. (2017b). Many useful relations for the spherical harmonic functions can also be found in, e.g., Dahlen and Tromp (1998). The spherical harmonics are a special case of the generalized spherical harmonics with Y m l = Y 0,m l in the convention that is used in this work.
The Lorentz force for a magnetic field of the form presented in Equation (12) can be written as where the radial functions of the radial component and the colatitudinal component are given by respectively. The coefficients are Ω N l = (l + N )(l − N + 1)/2 and γ l = (2l + 1)/4π. The last two factors of both equations are Wigner 3j symbols. Extensive lists of their properties can be found in, e.g., Edmonds (1960), Regge (1958) and Dahlen and Tromp (1998). We will be referring to only those properties of the Wigner 3j symbols that are listed in Appendix (E) of Kiefer et al. (2017b). Equations (15) and (16) are the vector spherical harmonic expansion coefficients of the Lorentz force for a toroidal magnetic field specified in Equation (12). They are derived in detail in Appendix B.
We expand aspherical perturbations to all stellar structural quantities in spherical harmonics (Dahlen and Tromp 1998;Lavely and Ritzwoller 1992;Woodhouse and Dahlen 1978): where the quantity Q can be the gravitational potential φ, density ρ, pressure p, or squared sound speed c 2 . We consider only zonal toroidal magnetic fields as perturbations. Hence, the azimuthal order of the spherical harmonic in the expansion (17) is set to 0. The ranges of the summation indices are determined by the configuration of the considered magnetic field: the indices s and s take the values of the model magnetic field. The index λ extends over all even values between 0 and s+s . This can be seen from properties of the Wigner 3j symbols (E29) and (E30c) in Kiefer et al. (2017b) and the fact that we only consider magnetic fields with even harmonic degree. Restricting the magnetic field to even harmonic degrees ensures antisymmetry about the equator which is generally observed for the Sun (Hathaway 2015). The aspherical perturbation to the gravitational potential δφ λ s,s (r) can be calculated by solving where G is the constant of gravitation, g 0 is the unperturbed gravitational acceleration, ρ 0 is unperturbed density, and R s,s ,λ and S s,s ,λ are the vector spherical harmonic coefficients of the radial and colatitudinal component of the Lorentz force. This equation is derived in Sweet (1950). Having obtained δφ λ s,s by numerically integrating Equation (18), we can calculate the perturbations in density, pressure, and squared sound speed: where Γ 1 is the first adiabatic exponent and c 2 0 is unperturbed sound speed. The derivatives of ln Γ 1 are supplied in the solar model that was used for this work (Christensen-Dalsgaard et al. 1996). We present brief derivations of Equations (19)-(21) in Appendix (D).
The general matrix element for the indirect effect is derived in Appendix (E) of Lavely and Ritzwoller (1992). We neglect terms due to rotation and ellipticity and transform the general matrix element from perturbations in bulk modulus κ 0 and density ρ 0 into perturbations in squared sound speed c 2 0 and density ρ 0 , see Equation (D40). This yields where K λ is the bulk modulus perturbation kernel and R (2) λ is the density perturbation kernel, which are listed in Appendix (C).

DISTURBING THE SUN
We modeled six different toroidal magnetic field distributions. Their radial profiles were modeled with Gaussians where B scale is a factor to scale the distribution to the desired maximum value. In Table 1 the locations of the maximum of the field distribution µ and the width in terms of the Gaussian standard deviation σ are given in the third and fourth columns. 1 The maximum values of the distribution are given in the fifth column. The sixth and seventh columns of Table 1 give the values of the plasma beta at the location of the maximum of the field distribution and at the photospheric level above the maximum of the distribution. Here, p is the gas pressure and B is the magnetic field strength. The plasma beta is the ratio of gas pressure to magnetic pressure. It must be noted that the values for β change with latitude as the spherical harmonics, with which the radial distribution are multiplied, see Equation (12), have a latitudinal dependence. This latitudinal dependence of the field distributions can be appreciated in the top panel of Figure 1, where magnetic field model A is shown in a meridional cut. The β values in Table 1 can thus be seen as minimal values.
Perturbations to the gravitational potential were neglected, i.e., we applied the Cowling approximation (Cowling 1941). This affects the sensitivity kernel R (1) λ (Equation (C28)) and the perturbations to the structural quantities (Equations (18)-(21)). 2 In our calculations, we used the standard solar model Model S by Christensen-Dalsgaard et al. (1996). Included in the version of the model we used is an extended set of variables, e.g., the derivatives ∂ ln Γ 1 ∂ ln p ρ and ∂ ln Γ 1 ∂ ln ρ p which are needed for the computation of the perturbation of the squared sound speed.
Model A is of harmonic degree s = 2, has its maximum at µ = 0.90 R with a width of σ = 0.04 R , and has a maximum field strength of B max = 50 kG, see the top panel of Figure 1. The maximum field strength is located at latitudes of ±45 • . The bottom panel of Figure 1 shows the resulting frequency shifts as a function of unperturbed mode frequency for modes with 4 ≤ l ≤ 148. To enhance clarity, we show only every fourth harmonic degree starting at l = 4. The frequency shifts are averaged over the azimuthal order m and are thus the mean shift of each multiplet. This shift is usually reported in studies of solar oscillation frequencies as a function of the level of activity (e.g., Broomhall 2017).
We calculated the general matrix elements for the direct and indirect effects separately. Thus, we can also examine the shifts caused by the two effects separately. In the top panel of Figure 2, the mean multiplet shifts of model A are shown for only the direct effect. In the bottom panel, the same is shown for the indirect effect. Note the different orders of magnitude of the shifts in the two panels.
Magnetic field model C, which is a strong field in the tachocline region with B max = 300 kG, is depicted in the top panel of Figure 3. We modeled this field with a harmonic degree of s = 2, at a depth of µ = 0.72 R , and with a width of σ = 0.04 R . The resulting frequency shifts as a function of unperturbed mode frequency for modes of harmonic degree 4 ≤ l ≤ 148 are shown in the bottom panel of Figure 3. Notice the small magnitude of the shifts compared to the field model A.
We chose the strength of this model because Antia et al. (2000) put an upper limit of 300 kG on magnetic fields in the tachocline region from analyses of mode splitting coefficients.
Model B is the same as model A, but with a lower maximum field strength. We will use this model in the next section to investigate differences between two configurations of the same geometry and depth but different strengths. Field configuration D has the same parameters for the radial function a(r) as models A and B, but has a harmonic degree of s = 4. Model E is a superposition of two fields, one with s = 2 and one with s = 4. Both fields are located at µ = 0.9 R with a width of σ = 0.04 R . The s = 2 component has a maximum field strength of B max = 50 kG, while the s = 4 component has a strength of B max = −30 kG. The cross-terms between the two configurations are included in the calculation of the general matrix elements. As can be seen in Figure 11 in Appendix E, due to  the superposition of the two fields, the maximum of the field is closer to the equator compared to model A.
Model F is a shallow field located at µ = 0.97 R with a width of σ = 0.01 R , with a maximum field strength of B max = 10 kG, and has harmonic degree s = 2. The multiplet shifts for this model are only of the order of some tens of nHz as can be seen in Figure 14 in Appendix E. The shifts caused by the indirect effect are much smaller still, being at the level of a few 10 −15 Hz. The shape of the frequency shifts as a function of mode frequency is rather different compared to the other models.   In Figure 4 the frequency shifts for models A and C are shown as functions of lower turning point of the modes. As can be seen in the bottom panel, only modes with lower turning point below or in the magnetized region around the tachocline experience a significant shift. Modes with turning points above the magnetic field are no longer disturbed by the direct effect. They can, however, still experience a weak shift due to the indirect effect, see lower row of panels in Figure 8 in Appendix E. As can be seen in the top panel of Figure 4, for model A, the maximum shifts are experienced by the  modes of the highest degree we calculated, as they have their lower turning point in the magnetized region.

DISCUSSION
The six magnetic field configurations we tested produce very distinct patterns in the frequency shifts as a function of unperturbed mode frequency. While models A, B, and D lead to a decrease of the multiplet frequency for a wide range of mode frequencies, the shift is generally positive for models C, E, and F. The multiplet shifts for field model A are presented in the bottom panel of Figure 1. The shifts for this model decrease for mode frequencies higher than ν ≈ 2 mHz and have their minimum at around ν ≈ 3.8 mHz. For even higher frequencies they increase and become positive for modes with ν 5 mHz. Ridges of modes with the same radial order are apparent. The pattern seen here is determined by the direct effect. As seen in Figure 2, the shifts caused by the indirect effect (bottom panel) are much smaller than those caused by the direct effect (top panel).
The mean multiplet shift caused by the indirect effect is negative for all models and all modes. The dependence of the shift as a function of frequency resembles the surface effect, see, e.g., Ball and Gizon (2014) or Basu (2016). For modes of low degree, as can be seen for the l = 4 modes, the lowest harmonic degree we considered here, the multiplet shift caused by the indirect effect is lower than for the next harmonic degrees we considered for models D and E. This peculiarity may help to rule out some field configurations if it is not observed for the Sun. We shall study the behavior of the low degree mode frequencies in a separate article as these are of special importance for asteroseismic studies of the magnetic activity of a star for which only modes up to l = 3 can be observed. We speculate that magnetic fields in the upper part of the solar convection zone can explain at least part of the discrepancy between model frequencies and observed frequencies.
The direct effect is the dominant contribution to the total shift for all model fields we investigated. The indirect effect is only important for more complex field configurations and when the magnetic field strength near the photosphere is so strong as to give a plasma beta β < 1. The shifts caused by the indirect effect are of the same order of magnitude as those caused by the direct effect for models D and E, as can be seen from Figures 10 and 12 in Appendix E. For the other models, which all have harmonic degree s = 2, the contribution from the indirect effect to the mean multiplet shift can well be neglected.
Generating plots for observed shifts as a function of the lower turning point of the modes, analogously to Figure 4, might give an indication of the location of magnetic field concentrations in the Sun. For this, the shifts of modes that have their lower turning point rather close to the surface have to be determined. If the shifts begin to decrease at a certain harmonic degree, as they do for the shifts for model C, their turning point can then be interpreted as the location of maximal magnetic field strength in this region. We note that such a location is not necessarily the location of the maximum field strength in the Sun, as deeper seated magnetic fields result in a smaller frequency shift even if they are stronger than more shallow magnetic fields. This can be seen by comparison of the top and bottom panels of Figure 4: model A, which has a maximum magnetic field strength of 50 kG located at 0.9 R , produces shifts at the µHz level. Whereas model C, which has a maximum magnetic field strength of six times that of model A, produces shifts that are only at the nHz level.
In Figure 5, the difference of the multiplet frequency shifts for model B and model A is shown as a function of unperturbed mode frequency. We show differences between the shifts of models B and A in the sense model B -model A, that is, the model with the weaker magnetic field minus the model with the stronger field. On the Sun, the frequency shifts are correlated with the level of magnetic activity. Hence, model B with a maximum field strength of 40 kG would correspond to the activity maximum and model A with B max = 50 kG would correspond to the activity minimum. In this picture, the toroidal magnetic field gets weaker in the ascending part of the cycle and reaches its minimum strength at the activity maximum. It then gets built up again by the solar dynamo  and reaches its maximum strength at the activity minimum. The large field strength then causes the activity to increase again as magnetic flux starts to rise to the surface.
In the top panel of Figure 6, we show the same as in Figure 5 but for a restricted frequency range of 1700-4000µHz. In the bottom panel of Figure 6, we show the results of Broomhall (2017, called B17 hereafter). She investigated the mean multiplet frequency shifts between the maximum of solar cycle 23 and the minimum of the same cycle from Global Oscillation Network Group (GONG) data. As can be seen, the differences between our two model magnetic field produce shifts in the multiplet frequencies that are strikingly similar to those reported by B17. Our results have the correct order of magnitude, while our frequency shifts are about 0.1 µHz higher than the solar shifts. This can, however, easily be corrected by adjusting the field strength of either model A or B. We find that the shifts of modes of low harmonic degree have a maximum at mode frequencies of ν ≈ 3800 µHz. As can be seen in Figure 5, they decrease for higher frequencies and even turn negative at about 5000 µHz. This behavior is not apparent in the result of B17. An extension of the measurement of the observed frequency shifts to mode frequencies above 4000 µHz may help to discern the model fields that best reproduce the shifts.
As can be seen in the last two columns of Table 1, the β values of the six models we consider here, cover a wide range. The tachocline field, model C, has the smallest β at its location of maximum magnetic field strength, but the highest value in the photosphere. Overall, the multiplet shifts for model C are very small, as can be seen in Figure 3. The shifts of model F, for which β(photosphere) = 154, are also only of the order of nHz. The other four models have β values at the photospheric level that are < 1. These models produce multiplet frequency shifts that are of the order of µHz. From this, we gather that the magnetic field needs to have an appreciable field strength in the near-surface layers if the mean multiplet shifts shall reach the µHz level as is observed for the Sun.
For the model field E, which is a superposition of two fields, we find that the shifts are positive for most modes, see Figure 12 in Appendix E. It is noteworthy that each of the two components of this model would lead to negative shifts on their own as can be seen for model A and model D. Also, the shifts are significantly different from zero already for mode frequencies below 2 mHz.

CONCLUSION
The search for the layers where strong magnetic fields are located in the interior of the Sun can be conducted with helioseismic forward calculations and inversions. In this article, we described our efforts to carry out forward calculations of the effect of toroidal magnetic fields on solar multiplet frequencies of acoustic oscillations. The theoretical framework we present here is also applicable to stellar models. We investigated six models for the magnetic field. A strong tachocline field, model C, produced only small shifts of the order of a few nHz. The maximum field strength of 300 kG is already at (Antia et al. 2000) or even high above the expected limit of magnetic fields in this region of the Sun (Arlt et al. 2007). With this, we can safely state that toroidal magnetic fields in the solar tachocline are not responsible for the observed frequency shifts over the solar cycle. Basu (1997) put an upper limit of 300 kG on the strength of toroidal magnetic fields concentrated below the convection zone base. This is the strength of magnetic field model C. Hence, the effect of fossil fields in the solar radiative zone with comparable strengths can also only be of the order of nHz.
Since we did not investigate poloidal magnetic field configurations, we cannot rule out poloidal fossil fields of this strength. We note again that the effect of near-surface magnetic field concentrations on acoustic oscillations are stronger by orders of magnitude and will thus likely impede the detection of fossil fields in the frequencies of solar acoustic oscillations. An investigation of the perturbation of eigenfunctions of low harmonic degree due to fossil fields is worthwhile and might put tighter constraints on the strength of such fields.
Shifts in multiplet frequencies of the Sun are measured between two differently magnetized states of the Sun. Commonly, these shifts are reported between an activity minimum and a following or preceding activity maximum. The difference of the multiplet shifts caused by our models A and B are found to be of the same order of magnitude as the observed shifts on the Sun. Also, the behavior of the shifts as a function of frequency and harmonic degree strongly resembles the solar values reported by B17. This might indicate that there are magnetic field distributions of the geometry and depth of our models A and B in the Sun. The location of the maximum field strength is at 0.9 R and they both have a field strength in the photosphere, which gives a plasma beta β that is smaller than unity. In general, we found that the β value in the very shallow layers must be of the order of unity to produce multiplet shifts that reach the observed µHz level.
A detailed study of splitting coefficients with the theory we presented here is expedient and will be conducted shortly. This will yield further insight into the depth, shape, and strength of the magnetic field concentration necessary to produce the observed frequency shifts and splitting coefficients caused by toroidal magnetic fields. It was shown by Schad et al. (2013) that analyses of perturbed eigenfunctions can be used to infer the solar meridional circulation. For this, the cross-coupling of modes has to be included, as self-coupling alone does not change the eigenfunction. Future integrated investigations of perturbed eigenfunctions, frequency shifts, and splitting coefficients, which can all be determined with the general matrix elements presented by Kiefer et al. (2017b) and in this article, are a promising tool to shed more light on the workings of the solar dynamo.
Acknowledgments -We are grateful to Jørgen Christensen-Dalsgaard for supplying the extended version of Solar Model S. We also thank Kolja Glogowski for computing the full set of solar oscillation eigenmodes as well as for helpful discussions, Vincent Böning for important comments and fruitful discussions, and our institute's internal referee Wolfgang Schmidt for reading the initial manuscript and his remarks, which helped to improve this article. We thank the anonymous referee for taking the time to review this paper. The research leading to these results received funding from the European Research Council under the European Unions Seventh Framework Program (FP/2007-2013)/ERC Grant Agreement no. 307117.

A. DECOMPOSITION OF A REAL VALUED VECTOR FIELD
In order to expand the Lorentz force in terms of spherical harmonics, we use that every real vector field u can be expanded in terms of vector spherical harmonics (Dahlen and Tromp 1998): The vector spherical harmonics, which are complete, are defined by where l is the harmonic degree, m is the azimuthal order with −m ≤ l ≤ m, Y m l is a spherical harmonic function, and we work in spherical geometry with the coordinates (r, θ, φ) radius, colatitude, and azimuth. The vector spherical harmonic coefficients in Equation (A1) are given by With Equations (A1)-(A7), the components of the vector field u are calculated as

B. EXPANSION OF THE LORENTZ FORCE
The Lorentz force, see Equation (11), for the toroidal configuration defined in Equation (12), is given by where a andâ are the radial profiles of the two magnetic field component, which are superposed. All dependencies were dropped. The Lorentz force (B11) will now be projected onto the vector spherical harmonics, as defined in Equations (A8)-(A10). Each vector component is treated separately.

B.1. The Radial Component
We concentrate on one magnetic field configuration with indices s, s , which is indicated by including the indices s, s to the notation of the vector spherical harmonic coefficients, e.g., R m s,s ,l .
Calculating the sum over different configurations in Equation (B22) then yields the total effect of the superposition.
With Equation (A8), the radial component of the Lorentz force can be written as where the vector spherical harmonic coefficient are We define the angular kernel where we used that Y 0,m l = Y m l , see Dahlen and Tromp (1998). The relations for the generalized spherical harmonics we use here can be found in Appendix (D) of Kiefer et al. (2017b). With Equation (D18) from Kiefer et al. (2017b), we can write The angular integral over the product of three generalized spherical harmonics can be calculated with help of Equation (C.198) of Dahlen and Tromp (1998). Making use of the properties of the Wigner 3j symbols (E30a)-(E30c) in Kiefer et al. (2017b), we find that µ = 0, as otherwise the angular integral vanishes. It can also be seen that the integrals over the first and last term in the bracket in Equation (B15) always vanish due to Equation (E30a) in Kiefer et al. (2017b). We thus find R s,s ,λ (r) = Ω s 0 Ω s 0 γ λ γ s γ s aâ r +â ∂a ∂r s s λ 0 0 0 where we made use of Equations (E27) and (E28) from Kiefer et al. (2017b) and we dropped the upper index on R as µ = 0 always holds.

B.2. The Azimuthal and Colatitudinal Components
Carrying out the same procedure for the azimuthal component as for the radial part leads to µ = 0. With Equations (A4), (A10), and (B11), we find that T 0 λ = 0 and hence F tor,φ = 0, as can be expected from Equation (B11). It remains to explicitly calculate the colatitudinal component of the Lorentz force: where the vector spherical harmonic coefficient is given by We define the two angular kernels The following properties and equations are then used to obtain the vector spherical harmonic coefficient: Equation (D18) from Kiefer et al. (2017b) to evaluate the colatitudinal derivatives; Equation (D16) from Kiefer et al. (2017b) with m = 0 to absorb the factor cot θ in the kernel A 2 ; Equation (C.198) from Dahlen and Tromp (1998) to evaluate the angular integral in Equation (B18); properties (E27)-(E30c) of the Wigner 3j symbols in Kiefer et al. (2017b). With all of this, we find where we dropped the upper index on S as µ = 0 always holds. The complete Lorentz force vector for a superposition of toroidal magnetic fields can thus be written as where the second summation extends over all even values between 0 and s + s due to properties of the Wigner 3j symbols (E28) and (E30c) from Kiefer et al. (2017b).

C. THE SENSITIVITY KERNELS
We reproduce the sensitivity kernels presented in Lavely and Ritzwoller (1992), Section (6c), which appear in the general matrix element for the indirect effect given in Equation (22) where K λ (r) is the bulk modulus perturbation kernel and R (2) λ (r) is the density perturbation kernel. The Woodhouse coefficients (Woodhouse 1980) are given by Applying the inverse chain rule and eliminating the radial derivative yields Equation (D34). The kernel for the indirect effect as given in Equation (90) of Lavely and Ritzwoller (1992) includes perturbations of the bulk modulus κ 0 and density ρ 0 . It can be transformed to perturbations in squared sound speed c 2 0 and density ρ 0 . The transformation is obtained by introducing perturbations in the defining equation of the squared sound speed c 2 0 = κ 0 /ρ 0 , subtracting the unperturbed equation, and linearizing in the perturbations. This yields δκ λ s,s (r) = ρ 0 (r)δc 2 λ s,s (r) + c 2 0 δρ λ s,s (r), which was used to obtain Equation (22) from Equation (90) in Lavely and Ritzwoller (1992). The aspherical perturbation to the bulk modulus can be expanded according to Equation (17). The perturbation in squared sound speed (Equation (D35)) is derived in Aerts et al. (2010), Section 3.6. Compared to Aerts et al. (2010), we neglect perturbations to the chemical abundances.