Red Giant Rotational Inversion Kernels Need Nonlinear Surface Corrections

Asteroseismology is our only means of measuring the rotations of stars in their interiors, rather than at their surfaces. Some techniques for measurements of this kind—“rotational inversions”—require the shapes of linear response kernels computed from reference stellar models to be representative of those in the stars they are intended to match. This is not the case in evolved stars exhibiting gravitoacoustic mixed modes: we show that the action of the asteroseismic surface term—systematic errors in the modeling of near-surface layers—changes the shapes of their inversion kernels. Corrections for the surface term are not ordinarily considered necessary for rotational inversions. We show how this may have caused previous estimates of red-giant envelope rotation rates from mixed-mode asteroseismic inversions to have been unintentionally contaminated by core rotation as a result, with errors comparable to the entire reported estimates. We derive a mitigation procedure for this hitherto unaccounted systematic error, and demonstrate its viability and effectiveness. We recommend this mitigation be applied when revising existing rotational inversions. Finally, we discuss both the prospects for applying such mitigation to the harder problem of inversions for stellar structure (rather than rotation), as well as the broader implications of this systematic error with regard to the longstanding problem of internal angular momentum transport.


INTRODUCTION
Almost all astronomical measurements of stellar properties are made only of their surfaces.Asteroseismology -the analysis and interpretation of stellar oscillations -is one of the very few ways we probe their interiors directly, through the use of so-called inversion techniques.In this paper, we focus on the rotational inverse problem (e.g.Schunker et al. 2016a,b), whose application to red giant solar-like oscillators (including but not limited to e.g Deheuvels et al. 2012Deheuvels et al. , 2014Deheuvels et al. , 2015;;Di Mauro et al. 2016;Fellay et al. 2021;Pijpers et al. 2021;Wilson et al. 2023) is a key lynchpin in our understanding of how angular momentum is internally transported and redistributed over the course of stellar evolution.
The nonrotating angular frequencies ω i of these oscillations, each associated with eigenfunctions ⃗ ξ i , are the normal-mode solutions to a Hermitian eigenvalue problem determined by the stellar structure (e.g.Unno et al. 1989;Gough 1993), expressed in terms of the wave operator L as L[ ⃗ ξ i ]+ω 2 i ρ ⃗ ξ i = 0.The response of these frequencies to small perturbations to the wave operator, L → L + ϵF , may be described by Rayleigh-Schrödinger perturbation theory through expansions of the form (e.g.Landau & Lifshitz 1965) under the convention that eigenfunctions are orthonormal with respect to the inner product: ⃗ ξ i , ⃗ ξ j = ⃗ ξ i • ⃗ ξ j ρ d 3 r = δ i j .For rotation rates Ω much lower than the critical frequency ω 0 ∝ GM/R 3 , individual modes of degree ℓ are split into multiplets of 2ℓ + 1 components each (e.g.Gough & Thompson 1990), with a multiplet splitting given to first order in Eq. (1) as Ong ignoring latitudinal differential rotation.Here K i is a sensitivity kernel (defined to be of unit integral, with the overall sensitivity instead specified by the normalisation factor β i ) for the i th mode.These integrals against sensitivity kernels in Eq. ( 2) are the diagonal matrix elements R ii of the angular momentum operator R (e.g.Lynden-Bell & Ostriker 1967), whose off-diagonal matrix elements also give rise to off-diagonal kernels by way of some quadratic form R acting on the mode eigenfunctions, defined in the sense that β i j K i j (r) Ω(r)dr = R i j = ⃗ ξ i , R ⃗ ξ j = R ⃗ ξ i , ⃗ ξ j Ω(r) dr.
(3) Such off-diagonal terms can be seen to enter Eq. (1) only at second and higher order.
In principle, this construction describes how the splittings of an actual star would relate to its (perfectly known) interior structure.However, the inversion kernels of an imperfect model may resemble those of the actual star, as long as these imperfections are small, in which case these integral expressions may also approximately apply to the kernels of such a reference model instead.The response of the mode frequencies to such small structural deviations may also take similar form to Eq. ( 2), integrating against sensitivity kernels only once over the stellar structure (e.g.Kosovichev 1999).This is a necessary condition for the notional "true" inversion kernels of the actual star to be well-described by those of the reference model used to approximate them.
One class of these rotational inversion techniques -that of optimal local averages (OLA: e.g.Backus & Gilbert 1968)operates on the principle that, given a sufficiently large number of observed modes, one may choose some linear combination of these sensitivity kernels, i c i β i K i , so as to suppress sensitivity except at a predetermined localisation radius r 0 (e.g.so that the summed kernel approximates a Gaussian of unit integral at that location).The corresponding linear combination of rotational splittings, i c i δω i , then yields an estimate of Ω(r 0 ), and thus permits Ω(r) per se to be inferred (by varying r 0 repeatedly; e.g.Christensen-Dalsgaard & Schou 1988).Operationally, then, these sensitivity kernels are found from stellar models constructed to closely match the observed properties of oscillating stars, with the inputs δω i to the inversion problem being the observed rotational multiplet widths.Similar principles underlie inversions for stellar structure (e.g.Basu et al. 2009), where the inputs are given by differences between the observed and model mode frequencies.For this purpose, the actual model mode frequencies are never directly used for this purpose -systematic error corrections are always applied, such as for the asteroseismic "surface term" (correcting modelling error in the near-surface layers).
Many variations on the OLA inversion procedure exist in the literature (e.g.Gough 1985;Pijpers & Thompson 1994;Rabello-Soares et al. 1999), with different ways of choosing coefficients {c i } (by optimising different merit/penalty functions) to yield summed kernels with various desirable properties.These all operate by penalising sensitivity of the summed kernel away from r 0 in different ways, while also, for example, minimising the expected observational error of i c i δω i , or maintaining a summed kernel of unit integral over the whole star.At their heart, they all nonetheless rely on the same fundamental assumption of linearity.This is also true of the other class of inversion techniques -the regularised least-squares method (RLS, e.g.Christensen-Dalsgaard et al. 1990), which fits constrained parameterisations of Ω(r) by forward-modelling multiplet splittings through Eq. ( 2) -also commonly applied in the literature.
While these linear inversion procedures require first-order accuracy for all modes used, the radius of convergence for perturbation theory varies from mode to mode.On this principle, Ong et al. (2021, hereafter OBR21) derive a nonlinearity threshhold for evolved solar-like oscillators like sub-and red giants.Should the surface term exceed it, such first-order perturbation analysis would then poorly describe their dipole modes of mixed gravitoacoustic character.This surface term acts as a further structural perturbation to a reference model, whose action may pass the pure p-modes through a forest of avoided crossings (e.g.their fig.11).Generically, such avoided crossings signify mode frequencies and eigenfunctions responding nonlinearly to structural or dynamical perturbations.Vanlaer et al. (2023) recently find that their presence in g-mode pulsators accompanies modifications to the shapes of the inversion kernels themselves, rendering structure inversions in those stars infeasible.This has historically not been a concern in helioseismology, in which astronomical context these techniques were first applied, nor for their application to main-sequence p-mode oscillators: the size of the surface term in those stars is small relative to the nonlinearity criterion of OBR21.In these stars, the surface term does not substantially modify the splitting widths given by Eq. (2), and thus is not thought to be relevant for rotational inversions.In any case, even where applied, existing prescriptions for surface-term corrections also operate only on mode frequencies, leaving eigenfunctions and inversion kernels unmodified.Current applications of rotational inversions to evolved sub-and red giants have inherited these helioseismic techniques with few modifications.While alternative, nonlinear, helioseismic inversion techniques exist (e.g.Marchenkov et al. 2000), these so far rely on specific features of analytic approximations to p-modes (Roxburgh & Vorontsov 1996), and are thus not yet suitable for use on mixed or g-modes.However, Li et al. (2023) report that typical sizes of the asteroseismic surface effect, as calibrated on Kepler red giants, are indeed large enough to advance modes along avoided crossings.Therefore, as this fundamental assumption of linearity no longer applies, the results of OBR21 immediately imply that the surface term changes the shape of inversion kernels in red giants too.Similar interference to Vanlaer et al. (2023) must thus afflict rotational inversions using these mixed modes.
We now examine how this may affect ongoing attempts at isolating seismic envelope rotation rates from red-giant dipole mixed modes.We moreover describe a well-defined mitigation procedure for the rotational inversion problem in particular.Finally, we conclude with some discussion regarding generalisations for the structural inversion problem, and broader implications.

CASE STUDY: ENVELOPE ROTATION OF KIC 9267654
To simplify the rotational inverse problem, red giants in particular are commonly treated in a two-zone model of differential rotation (e.g.Klion & Quataert 2017), with the core and envelope assumed to be rotating as more or less solid bodies (if at all).In general, fast rotation also induces nonlinear avoided crossings (e.g.Ouazzani et al. 2013;Deheuvels et al. 2017) when the off-diagonal matrix elements of Eq. ( 3) are large.However, these change the rotational splitting widths only to third order in perturbation theory (e.g.Ong et al. 2022), so we will, for the sake of argument, assume that the observed rotational splittings are small enough to remain well-described by the linear expression, Eq. ( 2).In that case, the rotational mutiplet width for a mode of mixed gravitoacoustic character in a red giant takes the form (Goupil et al. 2013) Here the angle brackets denote spatial averages over the gand p-mode cavities, corresponding roughly to the core and envelope -we leave the nature of this averaging underspecified -and ζ is a mixing fraction: ζ → 0 for a pure p-mode, and → 1 for a pure g-mode.In the asymptotic limit of high radial order, the normalisation constants take on limiting values of β p → 1, β g → 1 − 1/ℓ(ℓ + 1).
As an illustrative example, we examine a model of KIC 9267654: a rapidly rotating first-ascent red giant exhibiting such gravitoacoustic mixed modes at ℓ = 1, whose spectroscopic rotational broadening is greater than a putative envelope-averaged rotation rate implied by a combined seismic constraint from only ℓ = 2 and ℓ = 3 modes, suggesting rotational shear in the envelope (Tayar et al. 2022).This is a prima facie surprising result, as single-star evolutionary modelling instead favours close to solid-body envelope rotation.Given that these higher-degree modes are of significantly lower disc-integrated visibility than dipole modes, and penetrate less deeply into its interior, efforts are ongoing to measure envelope rotation rates from its dipole modes as well, to more strongly characterise this unusual rotational configuration.
Some of these ongoing studies of KIC 9267654 use the recently developed eMOLA technique (Ahlborn et al. 2022), which is in the OLA family of methods.Because of its simplicity and demonstrated robustness to noise in model tests, we examine it in detail here.However, our results concern the shapes of the kernels themselves, rather than the coefficients used to sum them, and so also impact all other inversion techniques (both OLA and RLS).Unlike standard OLA, eMOLA is intended only to isolate envelope rotation, rather than to localise rotation to a position r 0 , and derives inversion coefficients {c i } by optimising a penalty function constructed out of the cumulative sensitivity kernels β i (r) = r 0 β i K(r ′ )dr ′ .In addition to the usual contributions to the penalty function described above, the {c i } are additionally regularised so as to suppress i c i β i (r c ) at some radius r c situated just outside of the small radiative core.Such regularisation suppresses sensitivity of the summed kernel to the radiative core, owing to the highly oscillatory nature of the inversion kernels within the g-mode cavity rapidly cancelling each other out when summed and integrated.By comparison with Eqs. ( 2) and ( 4), this is equivalent to demanding that i c i β g,i ζ i be minimised.Mixed dipole modes in red giants specifically are of very high g-mode radial order, so the g-mode normalisation constants β g,i ∼ 1/2 do not vary significantly from mode to mode, and may be brought out of this sum as an overall multiplicative factor.
Since all the required properties of the inversion kernels used by the eMOLA procedure are effectively encapsulated in the mixing fraction ζ, we may recast the condition for optimality as minimizing the expected observational error of the inferred surface rotation rate, subject to constraints both that the summed kernel yields a normalised sensitivity of unity, and that the sum of mixing fractions is 0. In the language of Lagrange multipliers, we optimize an auxiliary penalty function where σ i are the observational errors on the rotational splittings.This gives a closed-form expression for the inversion coefficients c i as where a and b solve a 2 × 2 linear system of equations: We now apply these expressions to an illustrative mesa stellar model of KIC 9267654, which was the nearest simultaneous match to its effective temperature, metallicity, Gaia luminosity, radial and quadrupole p-mode frequencies (under the two-term surface-term correction of Ball & Gizon 2014), and dipole g-mode period spacing ∆Π 1 , from the model grid of Lindsay et al. (submitted to ApJ).Its dipole-mode frequencies and mixing fractions near ν max , computed using the pulsation code gyre (Townsend & Teitler 2013), are shown in blue in panels (a) and (b) of Fig. 1.We apply Eqs. ( 6) and ( 7) with values of σ i inversely proportional to a Gaussian envelope function situated at ν max , with width specified by the scaling relation of Mosser et al. (2012).This yields the linear combination of its rotational kernels shown with the blue curve in panel (c) of Fig. 1, which is insensitive to the core, as desired.Summed kernels constructed with the same coefficients, using rotational kernels for the reference (solid blue curve) and surface-perturbed (solid orange curve) models.The latter can be seen to retain sensitivity to the core.Repeating this exercise by cross-identifying modes between models with the lowest values of ζ can be seen to mitigate this only incompletely.(d): Residual senstivity to the core exhibited by the summed kernels of the surface-perturbed model, using coefficients of the reference model with nominal mode identification (orange curve), and under cross-identification in the neighbourhood of the most p-dominated modes (light gray curve).Without applying corrections, both can be seen to produce significant leakage.These can be seen to be suppressed by the use of coefficients obtained from the mitigation procedure where asymptotic mode-coupling calculations are performed after applying the surface correction to the underlying pure p-modes of the reference model (dark gray curve).The blue horizontal dashed line shows a notional value of 0. (Animation) The online animated version of this figure shows how the orange and gray objects of panels (a), (b), and (c) respond, as the stellar structure is linearly interpolated from the reference to the surface-perturbed stellar model over the course of 7.5 seconds.A vertical gray dotted line indicates the corresponding size of the surface perturbation in panel (d) during this process.
Let us now consider the corresponding summed kernel when the surface structure is slightly different, giving frequency differences from the reference model which are notionally those described by a surface-term correction.Since the true structure and rotational kernels of KIC 9267654 are not available for our inspection, we instead consider a stellar structure that is identical to that of our reference model, except for a small, artificially-induced perturbation to its near-surface layers.As in OBR21, we perturb the near-surface layers of the mesa model (bringing them out of hydrostatic equilibrium) as P(r) → P(r) × (1 + ∆(r)), with ∆(r) = A exp −(r/R − 1) 2 /2s 2 , and s = 0.002.Using a value of A ∼ 1 alters its pure p-modes by about 1 µHz at ν max , which is roughly in the range of values that Li et al. (2023) find to be needed for similar first-ascent red giants (cf.their fig.3).This is also roughly the size of the correction needed to bring the radial and quadrupole modes of this model into agreement with their observed values.We emphasise that in all other respects the inner structure of the perturbed model is identical to that of the reference; in particular the interior near-core layers, which determine the strength of coupling between the p-and g-modes, are untouched.
The resulting mode frequencies and mixing fractions (recomputed with gyre using identical settings as for the reference model) are shown in orange in panels (a) and (b) of Fig. 1.Importantly, as described in OBR21, the mixing fractions can be seen to change under the action of the surface term.This nonlinear interaction between the structure and mode frequencies also significantly changes the shapes of the inversion kernels compared to the reference model.Thus, if a linear combination of kernels of the reference model suppresses sensitivity to the core, that same linear combination of kernels from the surface-perturbed structure will fail to do so.This sum of the surface-perturbed kernels is depicted with the solid orange curve in panel (c) of Fig. 1.The corresponding linear combination of observed rotational splittings is thus nontrivially sensitive to core rotation, and ought not be interpreted as signifying only an isolated envelope rotation rate.In the depicted configuration, Fig. 1c shows that this adds up to 20% of the core rotation rate (on top of overall reduced sensitivity to the envelope of only 60%); we will refer to this as a 20% "core sensitivity leakage" in our following discussion.Given coefficients {c i,ref } of the reference model and mixing fractions ζ i,targ and sensitivity factors β i,targ of the target model, we compute this leakage L as A further complication to this situation is that, in actual observations, the identities of the modes (in particular their radial orders n pg ) are not known in advance.While we have assigned inversion coefficients from the reference model perfectly to modes with the same radial orders in the perturbed model (as is required by Eq. ( 2)), Fig. 1b shows that, under the action of a sufficiently large surface perturbation, the most p-dominated mixed mode in the reference model may not have the same radial order as that of the observed star.This being the case, one might argue that cross-identifying modes in assigning coefficients -so that the most p-dominated (and thus highest-amplitude) observed modes and rotational splittings are assigned the inversion coefficients associated with the most p-dominated modes in the reference model -might at least mitigate the issue.For the particular configuration shown, this is indeed the case.For illustration, we derive coefficients from the reference model for only the most pdominated mixed modes, and the 2 nearest ones on each side.
We then assign those coefficients to the most p-dominated mixed modes of the surface-perturbed model (indicated in gray in Fig. 1a and b) and likewise adjacent modes.The core sensitivity leakage of the resulting summed surface-perturbed kernel (gray curve of Fig. 1c) can be seen to be reduced, although not entirely eliminated, for this particular perturbed structure.
However, this introduces further problems, in the form of significant discontinuities and nonmonotonicity.We illustrate this with the orange and gray curves in Fig. 1d, which show how the core sensitivity leakage responds in both cases, as we smoothly interpolate the target model from the reference to the surface-perturbed structure.As the two structures diverge, the core sensitivity leakage under perfect mode identification increases in a monotonic fashion (orange curve), but that from cross-identifying modes does not, instead suffering unpredictable discontinuities with each discrete crossover of mode identification.The mitigation this affords is also in some cases actually counterproductive, leaking substantially more than direct mode identification.These are clearly undesirable properties of any mitigation procedure, and so such cross-identification is not a viable path forward.
3. MITIGATION STRATEGY So far, we have been restricted to qualitative analysis of an illustrative stellar model of a single star.However, this qualitative discussion indicates the existence of measurement systematics for envelope rotation, which we will now quantify.We consider the same set of evolutionary models as the hareand-hounds exercise of Ong & Gehan (2023), consisting of several evolutionary tracks of mesa models at different masses, solar composition, and solar-calibrated mixing length.Their p-and g-mode frequencies have been calculated using the prescription of Ong & Basu (2020).For each of these models, we compute the core leakages that would result from assigning the eMOLA coefficients of Eq. ( 6), as computed from the nominal mixed modes, to the modes and mixing fractions of the surface-corrected mixed modes, corrected according to the prescription of OBR21, and with coefficients for the pure p-modes of the two-term Ball & Gizon (2014) parameterisation computed from the calibration of Li et al. (2023).These are shown as a function of the relative g-mode density N = ∆ν/ν max 2 ∆Π 1 , a proxy for evolution up the red giant branch, in Fig. 2. We restrict our attention to 1 < N < 30, as the observed oscillations are primarily p-dominated modes at lower N, and mode mixing ceases to be an observational Ong concern for more evolved stars at higher N.While the amount of leakage appears to fluctuate significantly over the course of evolution (as the precise amount of leakage will depend on the relative positions of the underlying p-and g-modes), potentially even becoming negative, its severity can nonetheless be seen roughly to increase with evolution.This is in line with theoretical considerations, given that Li et al. (2023) finds the size of the surface term to be a roughly constant fraction of ∆ν, while the ratio of the nonlinearity threshhold of OBR21 to ∆ν decreases with increasing N. The cores of red giants are known to rotate more rapidly, by an order of magnitude or more, than their envelopes.Even a small, 10%, leakage of core sensitivity -well within the range of plausible values shown here -would then imply an amount of contamination comparable to the entire existing reported dipole-mode envelope rotation rates for red giants like KIC 9267654 (Triana et al. 2017), assuming these were not themselves contaminated in this fashion.Conversely, it is plausible that existing reported dipole-mode inferences of envelope rotation in such red giants -or even of counterrotation (e.g.Deheuvels et al. 2012Deheuvels et al. , 2014)), given the appearance of negative values for the core leakage in our modelling exercise -may be dominated by such contamination to begin with.Thus, the surface term is a hitherto unaccounted-for, but potentially dominant, source of systematic observational error in these rotational inferences.Clearly, there is a pressing need for mitigating measures.
We now outline a path forward.OBR21 described how standard surface-term corrections, applied to the pure p-modes of a stellar model (e.g. as computed through the π/γ decomposition scheme of of Ong & Basu 2020), yield surface-corrected mixed modes through mode-coupling calculations in which the mixed-mode eigenfunctions ⃗ ξ are expressed as linear combinations of pure p-and g-mode eigenfunctions.For mixed modes near each pure p-mode, only that p-mode contributes appreciably to this combination: to a good approximation, ⃗ ξ ∼ a p ⃗ ξ p + j a g, j ⃗ ξ g, j . (9) The coefficients of these linear combinations yield mixing fractions, as under unit normalisation ( ⃗ ξ, ⃗ ξ = 1), they satisfy Given that the rotational kernels are produced from a bilinear form (of the kind specified in Eq. ( 3)) acting on the eigenfunctions as βK = R[ ⃗ ξ, ⃗ ξ], we then have (10) Here we neglect cross-terms of the form R[ ⃗ ξ p , ⃗ ξ g ], as Ong et al. (2022) show that doing so, in this transformed basis of pure p-and g-modes, remains a good approximation even in the presence of nonlinear rotational avoided crossings in the natural basis of mixed modes.Angle brackets here represent a heuristic weighted average with respect to the squared mixing coefficients |a g, j | 2 .In words: we may well approximate the mixed-mode rotational kernels as linear combinations of pure p-and g-mode rotational kernels.Details for their computation, given a set of pure p-and g-mode eigenfunctions, are provided in Ong et al. (2022).We illustrate this construction in Fig. 3.
With this in hand, we may now propose our mitigation procedure.Given a reference model, one: 1. Computes the pure p-and g-mode frequencies and eigenfunctions from the reference model (e.g. using the prescription of Ong & Basu 2020), 2. Applies a surface term correction to the pure p-modes, 3. Performs mode-coupling calculations to obtain mixedmode frequencies and mixing fractions, and 4. Constructs mixed-mode rotational kernels of the form specified by Eq. ( 10).Subsequent inversion techniques are then to operate with respect to these reconstituted kernels.
Under the prescription of OBR21, step 3 requires solving for the matrix eigenvectors of a generalised Hermitian eigenvalue problem.However, the large rank of these matrices, particularly for evolved red giants, may render this procedure Mixed-mode rotational kernels may be constructed as linear combinations of pure p-and g-mode rotational kernels.Here, the cumulative integral for a mixed-mode rotational kernel of the surface-perturbed model (black dashed line), computed directly from its eigenfunction returned from gyre, is compared against a combination of that for an averaged pure g-mode rotational kernel (red solid line) scaled by the mixing fraction, with that for the nearest pure p-mode rotational kernel (gray solid line), scaled complementarily, returned from the reference model.prohibitively expensive.As an alternative, one may instead solve for the roots of the asymptotic characteristic function which arises from Jeffreys-Wentzel-Kramers-Brillouin (JWKB) analysis of the mixed-mode problem (Shibahashi 1979;Unno et al. 1989).Ong & Gehan (2023) derive a map from the characteristic function of the matrix problem to this JWKB characteristic function, as well as expressions for Θ p , Θ g , and q, in terms of the nonasymptotic pure p-and g-mode eigensystems.Given roots ν i of Eq. ( 11), the mixing fractions may then be found as the values of an asymptotic mixing function: The averaging over g-mode mixing coefficients a 2 g, j of Eq. ( 10) may also be replaced with an average over Brillouin-Wigner resonance factors, as a g, j (ν i ) ∝ 1/ ν 2 g, j − ν 2 i -we derive this explicitly in Appendix A. These approximate closed-form expressions, which hold best in the first-ascent-giant limit of low q, are much faster than the matrix calculation.
We show with the dark gray curve in Fig. 1d the core sensitivity leakage from summing kernels of the surface-perturbed model, using eMOLA coefficients obtained after applying our mitigation to the reference model, and using the above approximate JWKB expressions rather than solving the matrix eigenvalue problem.This simplified procedure can still be seen to efficaciously suppress such leakage without necessitating more expensive matrix calculations.We propose its adoption as a fast and easily implementable preprocessing step to the kernels used in existing code for rotational inversions with mixed modes.
4. DISCUSSION AND CONCLUSION Corrections for the surface term are not ordinarily considered necessary for rotational inversions.Even where such corrections have been applied, they have been thought to affect only the mode frequencies, and so all kernels that have been used for existing inversions remain those of the uncorrected modes in the reference models.We have now demonstrated how the action of the surface term changes the mixing fraction, and thus the shape of the rotational kernels, for mixed gravitoacoustic modes.As a result, this may have caused previous estimates of red-giant envelope rotation rates from mixed-mode asteroseismic inversions to have been unintentionally contaminated by core rotation.The potential scope of this issue is very broad, as it impacts effectively every existing mixed-mode rotational inversion.These unaccounted surface effects may explain the widely varying observational discrepancies -of comparable size to the reported values -between envelope rotation rates estimated from inversion techniques, compared to diagnostics derived from only the asymptotic theory without reference to stellar models (e.g.Goupil et al. 2013;Klion & Quataert 2017;Triana et al. 2017;Gehan et al. 2018;Mosser et al. 2018).While our above qualitative discussion was restricted to demonstrating core-rotation leakage into estimates of envelope rotation, it also conversely implies the potential existence of envelope-rotation leakage into existing estimates of core (differential) rotation.
We have described an explicit procedure to mitigate this issue.One first applies surface-term corrections to the mixed modes, via the generalisation proposed in OBR21, to obtain both surface-corrected mode frequencies, as well as surface-corrected mixing coefficients between the pure p-and g-modes.We have furthermore derived expressions for these mixing coefficients in terms of the commonly-used approximate JWKB formalism of mixed-mode coupling, to accelerate the matrix calculations otherwise required by the prescription of OBR21.One then uses these mixing coefficients to combine pure p-and g-mode rotational kernels, of the kind described in Ong et al. (2022), back into surface-corrected mixed-mode rotational kernels.
Similar considerations also apply to mixed-mode inversions for stellar structure, and may contribute to the so far inconclusive nature of attempts at this for stars more evolved than subgiants (Bellinger et al. 2021, and pers. comm.).However, the procedure that we have proposed here may not be directly applicable to the structure problem.Here, the shapes of the rotational kernels are changed only in response to perturbations to the stellar structure near the surface, which changes we reverse upon applying a nonlinear surface correction.Should the reference and target model differ in the stellar interior, and specificially in the evanescent region between the p-and g-mode cavities, then we may not necessarily be able to drop the p-g cross-terms that would appear in the structural equivalent to Eq. ( 10).Equivalently, a perturbation to these regions Ong may change the coupling strength between the two mode cavities, causing even the surface-corrected mixing fractions of the reference model to differ from that of the observed star.Better understanding the limitations of this procedure will be a crucial first step towards generalising it, and making further progress on the mixed-mode structure inversion problem.Corollarily, stellar models used for even the rotational inverse problem have to be constructed so that their gravitoacoustic coupling strengths match those of the observed stars.
While asteroseismic inversions are generally quite involved, our proposed mitigation requires only an additional preprocessing step to be applied to the inversion kernels -that is, only the penultimate step of a long analysis needs to be modified.This being the case, any revisions to these studies should be quick and unburdensome, so long as the original stellar models used for inversion remain available for reanalysis.To assist in these efforts, we have made our Python implementation of this mitigating prescription publicly accessible.This will in turn have immediate and astrophysically meaningful consequences.Seismic rotation measurements are how red giants inform the angular-momentum transport problem (e.g.Aerts et al. 2019;Fuller et al. 2019), and so such revisions may substantially modify our current understanding of it.We await these new findings with great excitement and anticipation.
We thank D. Huber, J. van Saders, and S. Basu, for constructive feedback and helpful discussions, and F. Ahlborn for insightful and enjoyable discussions about the eMOLA technique.JMJO acknowledges support from NASA through the NASA Hubble Fellowship grant HST-HF2-51517.001, awarded by STScI, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.Our Python implementation of the procedure described here, along with our input namelist files and the depicted reference mesa model, is freely available on Zenodo at doi:10.5281/zenodo.8433219.

APPENDIX
A. BRILLOUIN-WIGNER PERTURBATION THEORY Expressions like Eq. ( 1) emerge from standard Rayleigh-Schrödinger perturbation theory (e.g.Schrödinger 1926), which is used when both the perturbed eigenvalues and eigenvectors of a linear system are to be found simultaneously from the unperturbed eigensystem.However, when the perturbed eigenvalues are already known, one may instead make use of Brillouin-Wigner perturbation theory, which then yields expressions for the perturbed eigenvectors that are both more accurate (as the resulting perturbative expansion converges more rapidly), and also are far simpler, especially at higher order.We refer the reader to e.g.§3.5 of Killingbeck (1977) for a more thorough review regarding the properties of this approach to perturbation theory.However, since this is (to our knowledge) its first application in an asteroseismic context, we will also provide a rough overview of its main points here for quick reference.

(A2)
These expressions are understood to either be additionally scaled by an overall factor, or else to contain additional contributions parallel to |n (0) ⟩ at each succeeding order, if one wished to instead give unit normalisation for |n⟩.
In our application of Brillouin-Wigner perturbation theory, the perturbed frequency eigenvalues are the mixed-mode frequencies supplied from solving for the roots of the asymptotic characteristic function of Eq. ( 11).Ong & Gehan (2023) demonstrate that these correspond to the eigenvalues of the matrix eigenvalue problem, in the case of weak coupling between the pure p-and g-modes.Thus, for our purposes, our "unperturbed" system consists of the pure p-and g-modes, and the perturbation to the system, represented by V in Eq. (A2), is the coupling between them that produces mixed modes.This coupling has the further property of only being nonzero between p-and g-modes; that is, there is none between pure p-modes and other pure p-modes, or between pure g-modes and other pure g-modes.Finally, we may assume that the coupling strength may be treated as a smooth function of frequency.Thus, Eq. (A2) gives to leading order, and under unit normalisation, that for a g-dominated mixed mode, where the A i are normalisation factors such that j A i ω 2 p, j − ω 2 i 2 = 1.Conversely, for a p-dominated mixed mode, subject to the same kind of normalisation for the pure-g-mode component.It is this latter first-order expression which we use in Section 3.

Figure 1 .
Figure 1.In the presence of the surface term, a summed kernel constructed to suppress sensitivity to the core in the reference model will yield nontrivial sensitivity to the core with respect to a surface-perturbed set of modes.(a): Echelle diagram showing mixed-mode frequencies of a reference model (blue filled circles) of KIC 9267654, and one with a structural perturbation applied to its surface (orange open circles).Markers on the diagram are sized proportionally (by area) to the values of (1 − ζ) for the corresponding modes.Dashed lines indicate the locations of the underlying pure p-modes of the two models.The most p-dominated mixed modes of the surface-perturbed model are indicated with the filled gray circles.(b): Mixing fractions ζ for modes in the vicinity of the pure p-mode at ∼ 120 µHz.Coloured circles have the same meaning as in panel (a), and are joined with solid lines of matching colour.The most p-dominated surface-perturbed mixed mode is again indicated in gray.Dashed vertical lines denote the locations of the underlying pure p-modes of the models corresponding to each colour.(c):Summed kernels constructed with the same coefficients, using rotational kernels for the reference (solid blue curve) and surface-perturbed (solid orange curve) models.The latter can be seen to retain sensitivity to the core.Repeating this exercise by cross-identifying modes between models with the lowest values of ζ can be seen to mitigate this only incompletely.(d): Residual senstivity to the core exhibited by the summed kernels of the surface-perturbed model, using coefficients of the reference model with nominal mode identification (orange curve), and under cross-identification in the neighbourhood of the most p-dominated modes (light gray curve).Without applying corrections, both can be seen to produce significant leakage.These can be seen to be suppressed by the use of coefficients obtained from the mitigation procedure where asymptotic mode-coupling calculations are performed after applying the surface correction to the underlying pure p-modes of the reference model (dark gray curve).The blue horizontal dashed line shows a notional value of 0. (Animation) The online animated version of this figure shows how the orange and gray objects of panels (a), (b), and (c) respond, as the stellar structure is linearly interpolated from the reference to the surface-perturbed stellar model over the course of 7.5 seconds.A vertical gray dotted line indicates the corresponding size of the surface perturbation in panel (d) during this process.

Figure 2 .
Figure 2. eMOLA core sensitivity leakages computed for several tracks of solar-composition mesa models, assuming perfect mode identification, and assuming the surface term to take the parameterisation of Ball & Gizon (2014), with coefficients varying over evolution in the manner prescribed by Li et al. (2023).These are plotted as functions of the relative g-mode density N = ∆ν/ν max 2 ∆Π 1 : a proxy for evolution up the red giant branch.Models of each evolutionary track, indicated using filled circles coloured by mass, are shown connected with lines.Vertical positions of the markers indicate only the absolute value of the core leakage, with negative values indicated with additional open circles around individual markers.