1 Introduction↩︎

Magnetic fields are thought to play an important role in stellar structure and evolution, from inflating the radii of low-mass stars [1], [2], to affecting chemical transport by modifying convective overshoot [3], and transporting angular momentum via magnetic torques and magnetic instabilities (e.g., “Tayler–Spruit” dynamo, [4]–[8]). While many stars have observed surface magnetic fields [9], [10], much less is known about the magnetic fields deep within stellar interiors.

Asteroseismology [11] is one promising method to probe interior stellar magnetic fields. Weak magnetic fields produce asymmetric shifts in mode frequencies that can be calculated with perturbation theory [12]–[16], though [17] shows that perturbation theory can under-predict these shifts. Nevertheless, these types of asymmetric frequency shifts have recently been detected in tens of red-giant branch (RGB) stars [18]–[21].

The effects of strong magnetic fields on pulsation modes is still under debate. [22], [23] found that about 20 per cent of RGB stars have lower-than-expected dipole mode amplitudes, which [24] suggested was due to the presence of strong magnetic fields in their radiative cores. They argued that dipolar internal gravity waves (IGWs) would be scattered to higher wavenumbers by core magnetic fields, causing the gravity wave energy to be lost and depressing the mode amplitude. Others have used different techniques from perturbation theory to ray-tracing calculations to study the interaction of IGWs with core magnetic fields, finding that some of the wave energy is lost to magnetic waves, while some is reflected back up as IGWs, leading to a partial-suppression of the modes [25]–[28]. Indeed, [29] found that depressed dipolar modes still retain gravity-wave character, suggesting the IGW energy is not completely lost when they interact with magnetic fields. While no agreement has been reached regarding the efficiency of wave suppression due to strong core magnetic fields, wave suppression has been invoked to place limits on magnetic field strengths in a main-sequence star [30] and white dwarfs [31], as well as explain tidal dissipation in planetary systems [32].

To better understand the mechanism of IGW–magnetic field interaction, [33] ran two-dimensional (2D hereafter) numerical simulations in Cartesian geometry. They found that downward propagating IGWs perfectly convert to upward slow-magnetosonic (SM) waves, which develop smaller and smaller wavelengths as they propagate upward until they damp. The simulations were corroborated by a Wentzel–Kramers–Brillouin (WKB) analysis, in which the waves are assumed to be small-scale in the vertical (radial) direction, but large-scale in the horizontal (angular) direction. This contrasts with the ray-tracing analyses of [27], [28] which assume small horizontal scales, which may not be appropriate for studying dipolar IGWs. This work suggested that the gravity wave energy would be completely lost when interacting with a sufficiently strong magnetic field.

In the 2D wave problem there are two different wave modes: IGWs and SM waves. However, [34] pointed out that in three dimensions (3D hereafter) there exists a third wave mode, Alfvén waves (AWs), that may open new pathways for IGWs to interact with magnetic fields. They performed a similar WKB analysis as [33], but in spherical geometry and allowing for non-axisymmetric waves. They found that axisymmetric waves have the same properties as the 2D Cartesian waves described in [33]. For non-axisymmetric perturbations, in addition to the IGWs and SM waves, there also exist a continuous spectrum of resonant AWs. These resonant AWs were also described in [30], which modeled non-axisymmetric waves in spherical geometry using the same WKB approximation for the particular case of the main-sequence B star HD 43317. However, a major limitation of both of these previous works is that it is difficult to determine how this continuous spectrum affects wave propagation using WKB theory. This leaves open the question of whether or not interactions with the continuous spectrum could allow some IGW energy to reflect toward the stellar surface.

In this work, we revisit the problem of non-axisymmetric IGWs interacting with strong magnetic fields by solving the time-dependent MHD equations, and comparing the results to a WKB analysis. Following [33], we consider a simplified Cartesian model of the radiative zone of a RGB star. As in that previous work, we consider a magnetic field that varies both in \(x\) (representing latitude) and \(z\) (representing radius). While [33] assume no \(y\) dependence, here we assume the waves vary as \(\exp(i k_y y)\) with \(k_y\neq 0\). This allows for the presence of AWs, which can couple to the IGWs and SM waves. Despite the presence of this new wave mode, we find the same physical behavior as in the previous 2D simulations of [33]: there is perfect conversion from IGWs to magnetic waves (now a mix of SM waves and AWs), which damp when their wavenumbers become very large. Hence, we conclude that both axisymmetric and non-axisymmetric gravity wave energy is completely lost when interacting with a strong magnetic field.

The remainder of the paper is organized as follows. Section 2 details the Cartesian model and control parameters. Simulations conducted at two values of magnetic diffusivity are presented in Section 3; fine-scale features are interpreted using a simplified model of AW phase mixing. In Section 4, we characterize the large-scale wave modes by comparing our simulation results to an asymptotic solution found using a WKB approximation in the vertical direction. Finally, Section 5 summarizes our conclusions and discusses their astrophysical implications.

2 Problem setup↩︎

We consider a Cartesian model of the stably-stratified radiative core of a RGB star, with “latitude” \(x \in [0,L]\), “azimuth” \(y\), and “radius” \(z \in [0,L_z]\) with \(L_z = L/4\). Figure 1b shows a diagram of the domain over \(x\) and \(z\), roughly corresponding to the dashed region of the core radiative zone of a RGB star in Figure 1a.

The background state in our model is hydrostatic, with a stably-stratified density profile \(\rho_h(z)\), such that \(\partial_z \rho_h < 0\), and with a current-free magnetic field \(\boldsymbol{B}_0\), such that \(\boldsymbol{\nabla}\times\boldsymbol{B}_0 = \boldsymbol{0}\). Perturbations in velocity \(\boldsymbol{u} = u \boldsymbol{e}_x + v \boldsymbol{e}_y + w \boldsymbol{e}_z\), pressure \(p\), density \(\rho\), and the magnetic field \(\boldsymbol{b} = b_x \boldsymbol{e}_x + b_y \boldsymbol{e}_y + b_z \boldsymbol{e}_z\) are governed by the linearized magneto-Boussinesq equations [35]: \[\tag{1} \begin{equation}\tag{2} \overline{\rho} \partial_t \boldsymbol{u}+\boldsymbol{\nabla} p = -g \rho \boldsymbol{e}_z + \frac{1}{\mu_0}\left(\boldsymbol{\nabla}\times \boldsymbol{b}\right) \times \boldsymbol{B}_0, \end{equation} \begin{equation}\tag{3} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{equation} \begin{equation}\tag{4} \partial_t \rho = \frac{\overline{\rho} N^2}{g}w, \end{equation} \begin{equation}\tag{5} \partial_t \boldsymbol{b} = \boldsymbol{\nabla}\times(\boldsymbol{u}\times\boldsymbol{B}_0) + \eta \nabla^2 \boldsymbol{b}, \end{equation} \begin{equation}\tag{6} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{b}=0, \end{equation}\] where \(\overline{\rho}\) is the mean density (a constant), \(g\) is the gravitational acceleration, \(\mu_0\) is the magnetic permeability of free space, and \(\eta\) is the magnetic diffusivity. The dimensionless form of these equations used to derive our theoretical results are provided in Appendix 6.

We assume a constant Brunt-Väisälä frequency \(N = \sqrt{-g \overline{\rho}^{-1} \partial_z \rho_h}\). For simplicity, we ignore rotation, following [33] and [34]. This choice may be somewhat justified in the case of RGB stars, since the frequency of observed dipole modes \(\omega\) is much faster than the rotational frequency \(\Omega\) (\(\omega/\Omega \sim 20\), [23], [36]). In our simulations, magnetic diffusion is retained to regularize sharp features in the simulations associated with AWs (detailed in the next section), while viscous and radiative diffusion are neglected. Though this choice does not reflect the actual ordering of diffusivities in stars (typically, the radiative diffusivity is many orders of magnitude larger than \(\eta\)), it is a numerical convenience that allows the AWs to be resolved while minimally affecting the propagation of IGWs.

The background magnetic field is chosen to vary slowly in \(x\) and have vanishing vertical component at the “equator” (\(x/L = 0.25, 0.75\)): \[\boldsymbol{B}_0 = \mathcal{B} e^{-2\pi z/L}\left[\sin(2\pi x/L)\boldsymbol{e}_x + \cos(2\pi x/L)\boldsymbol{e}_z\right].\] The assumed exponential decay with \(z\) ensures that the field is current-free (i.e., \(\boldsymbol{J}_0 = \mu_0^{-1}(\boldsymbol{\nabla} \times \boldsymbol{B}_0) = \boldsymbol{0}\)) and thus stable to perturbations. Hence, \(\boldsymbol{B}_0\) may be thought of as the Cartesian equivalent of a dipole field. Field lines associated with \(\boldsymbol{B}_0\) are plotted in purple in Figure 1b.

An important control parameter for the waves is the Froude number, defined as \[\textit{Fr} = \frac{\text{wave inertia}}{\text{buoyancy}} =\frac{\omega}{N}.\] In the absence of magnetic fields, the Froude number controls the vertical scale of oscillations, as can be seen by rearranging the dispersion relation for pure IGWs: \[\label{eqn:igwdisprel} k_z^2 = (k_x^2 + k_y^2)(\textit{Fr}^{-2} - 1),\tag{7}\] where \(k_z\) is the “radial” wavenumber. In the limit of strong stratification, which we define as \(\textit{Fr} \ll 1\), the vertical wavelength scales as \(2\pi/k_z \sim l_z \equiv \textit{Fr} L\). This anisotropy is expected to hold for the suppressed dipole modes in RGB stars (\(\textit{Fr} \sim 10^{-1}\), [23], [37]) and is consistent with the Boussinesq approximation, for which vertical scales must be much smaller than the pressure scale height (the latter being comparable to \(L\)).

The presence of a magnetic field influences wave motion via the Lorentz force and Ohmic diffusion. The latter is characterized by the local-scale Lundquist number, \[\textit{Lu} = \frac{\text{wave inertia}}{\text{magnetic diffusion}} = \frac{\omega l_z^2 }{\eta} = \frac{\omega^3 L^2}{N^2 \eta},\] and is expected to have negligible effect on IGWs in the RGB core (\(\textit{Lu} \sim 10^{16}\text{–}10^{20}\), [23], [37], [38]).

The strength of the magnetic field is measured in terms of the dimensionless “magneto-gravity ratio” \[\Gamma = \frac{\text{Lorentz force}}{\text{wave inertia}} = \left(\frac{\omega_{\textit{MG}}}{\omega}\right)^2,\] where \[\omega_{\textit{MG}} = \left( \frac{\mathcal{B} N}{L\sqrt{\mu_0 \overline{\rho}}}\right)^{1/2}\] is the magneto-gravity frequency [24]. For the stellar cores of interest, we cannot determine the value of \(\Gamma\) a priori since \(\mathcal{B}\) is not known at the outset. Instead, we follow [24] in assuming that the background magnetic field is strong enough to influence IGWs but not so strong as to inhibit horizontal (\(x\),\(y\)) perturbations, \(\boldsymbol{u}_H\), \(\boldsymbol{b}_H\). Accordingly, we assume that \[\lvert \partial_t \boldsymbol{u}_H \rvert \gtrsim \left\lvert \frac{\boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{e}_z}{\mu_0 \overline{\rho}} \partial_z \boldsymbol{b}_H\right\rvert.\] Using \(l_z = \omega L/N\) as the scale for the vertical gradient, \(\mathcal{U}\) as the (arbitrary) scale of \(\boldsymbol{u}_H\), and \(\mathcal{U} \sqrt{\mu_0 \overline{\rho}}\) as the scale of \(\boldsymbol{b}_H\), this implies that \[\omega \mathcal{U} \gtrsim \frac{\mathcal{B}}{\mu_0 \overline{\rho}} \frac{1}{l_z} \mathcal{U} \sqrt{\mu_0 \overline{\rho}},\] and therefore \[\Gamma = \frac{\mathcal{B} N}{\omega^2 L \sqrt{\mu_0 \overline{\rho}}} \lesssim 1.\]

The magneto-gravity ratio \(\Gamma\) used in this work is similar to the “depth parameter” \(a\) in [34]. However, \(a\) is proportional to the local magnetic field strength (and may be used as a vertical coordinate) whereas \(\Gamma\) is proportional to the strength \(\mathcal{B}\) of the background field at the base of the domain (and is thus a constant). In our problem, this is enough to characterize \(\lvert \boldsymbol{B}_0 \rvert\) everywhere since \(\lvert \boldsymbol{B}_0 \rvert = \mathcal{B} \exp(-2\pi z/L)\) .

To study the interaction of IGWs with magnetic fields, we numerically solve initial value problems (IVPs) at \(\textit{Fr} = 0.025\) and \(\Gamma = 0.1\), in which waves are forced at a fixed frequency \(\omega\) and horizontal wavenumbers \(k_x = k_y = 2 \pi/L\). Though we cannot reach astrophysical values of \(\textit{Lu}\), we show that the relevant physics already occur at more moderate values of \(\textit{Lu}\) by comparing two cases with \(\textit{Lu} = 6.25 \times 10^4\) (IVP I) and \(\textit{Lu} = 6.25 \times 10^5\) (IVP II).

We assume that all dependent quantities can be expressed as \(q(x,y,z,t) = \Re\{q'(x,z,t) \exp({i k_y y})\}\) (so \(\partial_y \equiv i k_y\)), and we solve the remaining partial differential equations in \(x,z,t\) using the Dedalus pseudospectral code [39] with (\(N_x\), \(N_z\)) complex Fourier modes in both directions. A two stage, second-order implicit–explicit Runge–Kutta scheme [40] with a uniform timestep \(\Delta t\) is used to evolve the flow from an initially quiescent state. The timestep, resolution, and other simulation parameters are reported in Table [tab:simparams].

To excite “dipolar” IGWs with fixed frequency \(\omega\) and latitudinal wavenumber \(k_x = 2 \pi/L\), we add a forcing term \[F = \exp\left[-\frac{(z-z_0)^2}{\Delta z^2}\right]\exp\left[i(k_x x + k_y y - \omega t)\right],\] to the right-hand side of the density equation (4 ), where \(\Delta z = L/240\) is the vertical width of the forcing and \(z_0 = 9 L/40\) is the forcing height (dotted line in Figure 1b). We set the azimuthal wavenumber of the forcing to \(k_y=k_x=2\pi/L\) by analogy with non-axisymmetric dipole modes, which have \(\lvert m \rvert = \ell = 1\).

In order to avoid unphysical reflections of the waves at the top and bottom boundaries of the domain, we add damping layers near \(z = 0\) and \(z = L_z\) (hatched regions in Figure 1b). These layers are implemented by adding linear damping terms \(-D \boldsymbol{u}\), \(-D \rho\), and \(-D \boldsymbol{b}\) to the right hand sides of (2 ), (4 ), and (5 ), respectively. The damping rate \[D(z) = \frac{\omega}{2}\left[ \tanh\left(\frac{z - L_z + s}{\Delta z}\right) + \tanh\left(\frac{s - z}{\Delta z}\right) + 2\right]\] ensures that perturbations are damped on timescales longer than \(\omega^{-1}\). The size of the layers is given by \(s = L/80\).

The use of Fourier modes in \(z\) requires the background magnetic field to be periodic. We achieve this by multiplying \(\boldsymbol{B}_0\) by \(-\frac{1}{2}[ \tanh\left(({z - L_z + s_B})/{\Delta z_B}\right)\) \(+ \tanh\left(({s_B - z})/{\Delta z_B}\right)]\), where \(s_B = s/3\) and \(\Delta z_B = \Delta z/4\) so that the field is tapered to \(\boldsymbol{0}\) within the damping layers at the top and bottom of the domain.

The simulations are integrated until \(t = 1000 \omega^{-1}\), long after the time at which the mean wave energy per unit volume \[\label{eqn:meanenergy} E = \frac{k_y}{2\pi L_x L_z}\int_0^{L_z}\int_{-\pi/{k_y}}^{\pi/{k_y}}\int_0^{L_x} \mathcal{E}\mathrm{d}x\mathrm{d}y\mathrm{d}z,\tag{8}\] saturates (\(E\) reaches 99 per cent of its final value at least 75 wave periods before the end of each simulation). In (8 ), the local wave energy density is \[\mathcal{E} = \frac{1}{2}\left(\overline{\rho}{ \boldsymbol{u}\cdot\boldsymbol{u}} + \frac{g^2 \rho^2}{\overline{\rho} N^2}+\frac{\boldsymbol{b} \cdot \boldsymbol{b}}{\mu_0} \right).\]

3 Simulation results↩︎

Figure 2 shows snapshots of IVP I (\(\textit{Lu} = 6.25 \times 10^4\)) in the equilibrated state. Panel a presents the latitudinal flow \(u\) in the \((x,z)\) plane located at \(y=0\); the corresponding azimuthal flow \(v\) is shown in panel d. The domain is vertically stretched to better show the rapid oscillations in \(z\), as anticipated from dimensional analysis (\(l_z = \textit{Fr} L\), \(\textit{Fr} = 0.025 \ll 1\)). Hatched regions indicate the damping layers.

Since the \(x\) and \(z\) components of the background magnetic field \(\boldsymbol{B}_0\) have distinct parities, the sine- and cosine-parity components of \(u\) and \(v\) behave differently upon interacting with \(\boldsymbol{B}_0\); it is therefore instructive to analyze these components separately. To satisfy incompressibility, \(u\) and \(v\) must have opposite parity in \(x\): the sine component of \(u\) (shown in panel b) is associated with the cosine component of \(v\) (shown in panel e), and the cosine component of \(u\) (panel c) is associated with the sine component of \(v\) (panel f). Hereafter, we shall refer to waves associated with sine-parity (cosine-parity) latitudinal flow \(u\) as sine-parity (cosine-parity) waves.

The forcing (dotted line at \(z_0/L = 0.225\)) excites IGWs with both sine and cosine parity (checkerboard features in Figure 2b,c,e,f). At early times in the simulation, these IGWs travel downwards into a region of increasingly strong magnetic field. Below \(z/L \approx 0.05\) (still many wavelengths above the bottom damping layer), the field becomes too strong to permit wave propagation, and now causes the down-going waves to become evanescent instead (see Section 4 for more detail). Above \(z/L \approx 0.05\), the strong magnetic field causes the waves to reflect and propagate back in the upwards direction. These up-going waves (of both sine and cosine parities) have different horizontal structures than the IGWs forced at the top of the domain.

The latitudinal flow \(u\) associated with the sine-parity up-going wave (Figure 2b) develops shorter and shorter vertical scales as the magnetic field weakens with height, until it vanishes near \(z/L \approx 0.17\). This apparent cut-off height resembles the SM cutoff heights in the axisymmetric case [33], as discussed further in Section 4. However, for \(k_y \neq 0\), the coupling of the latitudinal flow \(u\) to the azimuthal flow \(v\) introduces additional complexity. The sine-parity up-going wave is associated with azimuthal flow \(v\) (shown in panel e) that exhibits large-scale ridges centered at \(x/L = 0, 0.5, 1\). The amplitude of these ridges decays gradually above the cutoff height near \(z/L \approx 0.17\) (in contrast to the concurrent sharp decay in \(u\)). At the cutoff height, energy in the up-going wave is transferred from \(u\) to \(v\), coinciding with the growth of fine-scale tilted features in \(v\) on either side of each ridge.

The latitudinal flow \(u\) associated with the cosine-parity up-going wave (Figure 2c) also vanishes at a cutoff height located this time around \(z/L \approx 0.1\). At this point, energy is transferred to the sine-parity azimuthal flow (shown in panel f), which develops features with extremely fine horizontal scales.

The axisymmetric simulations of [33] also contain small-scale oscillations, for which there was no physical interpretation in that work. Those features are extremely limited in latitudinal extent (concentrated around \(x/L = 0.25, 0.75\) where the vertical component of the background field vanishes) in contrast with the fine-scale features in Figure 2e,f. As discussed further below, we interpret these features (in the present simulations) as AWs.

AWs are known to develop small horizontal scales as a consequence of phase mixing (due to the variation of Alfvén speed with \(x\)) and rapidly damp via Ohmic diffusion [41], [42]. To show how phase mixing occurs in our simulations, we present a ‘toy’ model (after [42]), involving a purely vertical background magnetic field \(B_{0 z}(x)\boldsymbol{e}_z\) that varies with \(x\). We consider open field lines and an axisymmetric forcing that displaces fluid at the base of the domain by \(\xi_0(x)\exp(-i\omega t)\boldsymbol{e}_y\), where \(\xi_0(x)\) is an arbitrary displacement profile.

In the absence of diffusion, each magnetic field line can move independently in the \(y\) direction and thus supports AWs associated with azimuthal fluid displacement \[\label{eqn:toydisplacement} \xi_y = \frac{v}{-i\omega} = \frac{b_y}{i k_A B_{0 z}} = \xi_0(x)\exp[i k_A(x) z-i\omega t]\tag{9}\] where \[\label{eqn:toyalfvenwavenumber} k_A(x) = \frac{\omega \sqrt{\mu_0 \overline{\rho}}}{ \lvert B_{0 z}(x) \rvert}\tag{10}\] is the local Alfvén wavenumber.

To illustrate the resultant phase mixing process, we consider a basal displacement profile \(\xi_0(x) = \sin(2\pi x)\) and background magnetic field with \(B_{0 z}^2 = \cos^2(2\pi x) + 0.01\) to ensure that \(k_A(x)\) remains bounded for all \(x\). Figure 3 shows magnetic field lines \(B_{0 z}\boldsymbol{e}_z + b_y\boldsymbol{e}_y\) (purple curves) associated with the forced AWs. Profiles of the azimuthal displacement \(\xi_y\) over \(x\) are overlaid in orange at four different heights. Since stiffer field lines (higher \(\lvert B_{0 z}(x) \rvert\), darker purple curves) oscillate with longer wavelengths (lower \(k_A\)), field lines that are in phase at one height become out of phase higher up. This process of phase mixing produces small horizontal scales in \(\xi_y\) (and thus \(v\)) that become finer with height (compare the top and bottom orange curves). The finest scales are concentrated where \(k_A\) changes most rapidly with \(x\), around minima in \(\lvert B_{0 z}(x) \rvert\) (at \(x = 0.25, 0.75\)).

Phase mixing in IVP I produces similar features in the azimuthal flow. For the cosine-parity wave, Figure 2f shows fine scales in \(v\) above \(z/L \approx 0.15\), concentrated around \(x/L = 0.25, 0.75\) (near which \(\lvert \boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{e}_z \rvert\) changes rapidly), as in the toy problem. However, our simulations contain further complexity: the background field is not purely vertical, magnetic diffusion is present, and the waves are non-axisymmetric such that azimuthal and latitudinal motions are coupled via continuity. Further, sharp features in \(v\) (not due to phase mixing) are already present in the up-going large-scale wave modes below their cutoff heights.

For the cosine-parity wave, these sharp features (the tilted “shingles’’ in Figure 2f below \(z/L \approx 0.1\) and away from \(x/L = 0.25, 0.75\)) are located at critical layers where the large-scale wave mode resonates with the Alfvén spectrum, as will be demonstrated in Section 4. At the cutoff height (near \(z/L \approx 0.1\)), the large-scale wave mode acts analogously to the basal forcing in the toy problem, transferring energy to the phase-mixing AWs. In the toy problem, the sinusoidal structure of the forcing is imparted to the AWs. Similarly, in the simulation, the AWs just above \(z/L \approx 0.1\) retain the critical layer shingles present in the large-scale wave mode.

The large-scale up-going wave with sine parity also converts to AWs above its cutoff height (near \(z/L \approx 0.17\)). As in the cosine-parity case, features of the large-scale wave (i.e., the ridges centered at \(x/L = 0,0.5,1\) in Figure 2e) persist in the azimuthal flow associated with the AWs and gradually decay with height. The effect of phase mixing is weaker than before since \(v\) vanishes at \(x/L = 0.25, 0.75\) for the sine-parity waves. Nonetheless, thin tilted shingles associated with phase mixing are still visible next to the ridge features in Figure 2e, and they grow finer towards \(x/L = 0.25, 0.75\) as expected.

Fine-scale AWs of both parities are efficiently damped by magnetic diffusion. In order to isolate the effects of resistivity on the AWs and large-scale modes of oscillation, we run a second IVP (IVP II) and decrease the magnetic diffusivity so that \(\textit{Lu} = 6.25 \times 10^5\).

Figure 4 shows \(u\) and \(v\) decomposed by parity for IVP II. We see that the large scale features observed in IVP I are still present in IVP II, and are mostly unaffected by the decrease in magnetic diffusivity. Similarly, AWs emerge as before at \(z/L \approx 0.17\) and \(z/L \approx 0.1\) for the sine and cosine parities (Figure 4e and f), respectively. However, \(v\) develops significantly finer scales in IVP II than in IVP I, since diffusion sets the minimum length scale of AWs [43].

In the next section, we characterize the large scale features as wave modes found via a WKB approximation in the vertical direction. The following WKB analysis elucidates the mechanism of wave conversion and explains why the large scale modes have the same behavior in IVP I and II.

4 WKB analysis↩︎

In this Section, we use a WKB analysis to derive an expression for the large-scale oscillation modes in the limit of small \(\textit{Fr}\). The approximation takes advantage of the separation between the short vertical scale of the oscillations (\(O(\textit{Fr} L)\)) and the slow variation of the background magnetic field (\(O(L)\)). We approximate each field as, e.g., \[\label{eqn:waveansatz} p \sim A(z) p_0(x,z)\exp\left(i\theta(z) + i k_y y - i \omega t\right),\tag{11}\] where \(k_y = 2\pi/L\) and \(\omega\) are fixed by the forcing. The local scale of oscillations is described by the vertical wavenumber \(k_z(z) = \partial_z \theta\). Crucially, the WKB approximation assumes that \(\theta(z)\) varies rapidly with \(z\) over \(O(\textit{Fr} L)\) scales while \(k_z(z)\), \(p_0(x,z)\), and \(A(z)\) vary slowly with \(z\) on the system scale \(L\). The amplitude function \(A(z)\) is common to all the fields (\(\boldsymbol{u}, p, \rho, \boldsymbol{b}\)), and is found by imposing a solvability condition on the next-to-leading order equations. The details of the WKB procedure and amplitude equation (29 ) may be found in Appendix 7.

We assume that the strength of the magnetic field is commensurate with the degree of stratification such that \(\Gamma\) remains \(O(1)\) as \(\textit{Fr} \to 0\). In this distinguished limit, we show in Appendix 7 that the WKB approximation reduces the linear equations (1 ) to a one-dimensional generalized eigenvalue problem in \(x\) that depends parametrically on \(z\) through \(\boldsymbol{B}_0(x,z)\). Ignoring the effects of magnetic diffusion, the problem can be reduced to a single equation for the pressure eigenfunction \(p_0(x,z)\) and eigenvalue \(k_z^2(z)\): \[\label{eqn:idealpressure} \frac{k_z^2}{N^2} p_0 + \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot}\left(\frac{ \hat{\boldsymbol{\nabla}} p_0}{\omega^2 - v_{A z}^2 k_z^2}\right) = 0,\tag{12}\] where \(\hat{\boldsymbol{\nabla}} = \boldsymbol{e}_x \partial_x + \boldsymbol{e}_y i k_y\) is the transformed horizontal gradient operator and \[v_{A z}(x,z) = \frac{\boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{e}_z }{\sqrt{\mu_0 \overline{\rho}}}\] is the local Alfvén speed. (Note that this eigenproblem can be solved independently at each height \(z\) in the domain). The pressure equation (12 ) is the Cartesian analog of the spherical eigenvalue equation in [34].

If \(v_{A z} \equiv 0\), then (12 ) reduces to the (pure) IGW dispersion relation in the \(\textit{Fr} \to 0\) limit (cf. Equation 7 ): \[\label{eqn:lowFrigwdisprel} \frac{k_z^2}{N^2} - \frac{k_x^2 + k_y^2}{\omega^2} = 0,\tag{13}\] for which the eigenfunctions are Fourier modes \(p_0 = \exp(i k_x x)\) and the eigenvalues (\(k_z^2\)) are independent of \(z\). In our model, \(v_{A z} \neq 0\) and \(\partial_z v_{A z} \neq 0\), so both the eigenfunctions \(p_0(x,z)\) and the vertical wavenumber \(k_z(z)\) must vary (slowly) with \(z\).

4.1 Axisymmetric modes↩︎

Figure 5a shows the real part of the vertical wavenumber \(k_z\) versus height \(z\) for the axisymmetric (\(k_y = 0\)) case (similar to the constant \(N\) case analyzed by [33]). The eigenvalues \(k_z^2\) are found as described in Appendix 8.1. Only the wavenumbers of “dipolar” modes (i.e., those for which \(\lvert p_0\rvert\) has two zeros on \(x \in [0,L)\) at each \(z\)) are included. Arrows are oriented along the wavenumber branches according to the sign of the WKB (vertical) group velocity, which is computed using (30 ).

Figure 5a shows 2 overlapping branches of \(k_z\) corresponding to the sine-parity IGW (IGW-0, yellow curve) and the cosine-parity IGW (IGW-1, orange curve), respectively. As before, we refer to the modes by the parity of the associated latitudinal flow, \(u_0\). Both IGW-0 and IGW-1 branches closely follow the prediction for pure IGWs (\(k_z = \textit{Fr}^{-1}\sqrt{k_x^2 + k_y^2}\), vertical gray dash-dotted line) near the top of the domain, but refract to higher wavenumber as they propagate down into regions of greater magnetic field strength.

At \(z/L \approx 0.0125\), IGW-1 reaches a turning point, below which it leaves an an evanescent tail (evan.-1, dotted orange curve in Figure 5a) with \(\Im\{k_z\} < 0\). Note that for the given magnetic field strength (\(\Gamma = 0.1\)), the corresponding turning point for IGW-0 is outside the domain at \(z/L = \ln(2\pi \Gamma/1.08)/(2\pi) \approx -0.0862\) [33].

Above the IGW-1 turning point, there is a connecting wavenumber branch corresponding to an up-going cosine-parity SM mode (SM-1, purple curve in panel a). Since the restoring mechanism for SM waves is dominated by the Lorentz force rather than buoyancy (as for IGWs), the properties of SM-1 are intrinsically tied to the background magnetic field. The wave refracts again as it propagates up to regions of lower magnetic field strength. The upward path of SM-1 is bounded by a magnetic cutoff height (Figure 5a, horizontal gray dashed line) at \(z/L = \ln(6\pi \Gamma)/(2\pi) \approx 0.101\) [33], near which \(k_z\) becomes infinitely large and the wave damps for arbitrarily small diffusivity.

[34] found an important difference between these behaviors (in a Cartesian domain) and the corresponding behavior of IGWs and SM waves in spherical geometry with a dipolar background field: for axisymmetric dipole modes (\(\ell=1\), \(m=0\)), they find that the IGW approaches a magnetic cutoff height from above and never reaches a turning point. However, they show that \(\ell > 1\), \(m=0\) modes exhibit similar IGW\(\rightarrow\)SM turning points as in the Cartesian analysis of [33] and Figure 5 (a). Moreover, we expect the inconsistent behavior for the \(\ell=1\), \(m=0\) case to have little bearing on non-axisymmetric dipole modes since [34] also find IGW\(\rightarrow\)SM turning points for \(\ell = |m|=1\), as discussed in Section 4.2.

As [34] further note, the ideal magnetogravity eigenproblem (12 ) has an internal singularity at \[\label{eqn:AWdisp} \omega^2 - v_{A z}^2 k_z^2 = 0,\tag{14}\] which is also present in the polarization relation for \(v_0\) (the azimuthal velocity eigenfunction): \[\label{eqn:vpolrel} v_0 = \frac{k_y \omega p_0}{\overline{\rho}(\omega^2 - v_{A z}^2 k_z^2)}.\tag{15}\] Since the local Alfvén speed \(v_{A z}\) varies continuously with \(x\), this singularity is associated with a continuous spectrum of AW modes (gray shaded region in Figure 5a) with \[\label{eqn:alfvenbdry} \lvert k_z \rvert \geq k_{A B}(z) = \frac{\omega}{\max_x \lvert v_{A z}(x,z)\rvert} =\frac{1}{\Gamma l_z} \exp\left(\frac{2\pi z}{L}\right)\tag{16}\] satisfying (14 ) somewhere in the domain. The separation between the IGW branches (\(k_z \sim l_z^{-1}\)) and the Alfvén continuum (\(k_z \sim \Gamma^{-1} l_z^{-1}\)) decreases inversely with the magneto-gravity ratio \(\Gamma\).

When \(k_y = 0\), the singularity in (12 ) may be removed after integrating the equation with respect to \(x\). Then, (12 ) and (15 ) yield decoupled dispersion relations involving the meridional (\(x\), \(z\)) streamfunction \(\psi_0 = (\omega k_z\int_0^x p_0 \mathrm{d}x')/(\overline{\rho} N^2)\) and the azimuthal flow \(v_0\), respectively: \[\begin{equation} \omega^2 k_z^2 \psi_0 - v_{A z}^2 k_z^4 \psi_0 + N^2 \partial_x^2 \psi_0 = 0, \end{equation} \begin{equation} (\omega^2 - v_{A z}^2 k_z^2)v_0 = 0. \end{equation}\] In this case, each AW eigenfunction \(v_0\) may be expressed as the sum of Dirac delta functions \[\label{eqn:AWaxisymmode} v_0 = \frac{1}{N_c}\sum_{n=1}^{N_c} \delta(x-x_{c,n})\tag{17}\] centered at critical latitudes \(x_c\) where \(\omega^2 - v_{A z}^2 k_z^2 = 0\). Physically, (17 ) and (16 ) mean that each magnetic field line oscillates in the azimuthal direction with a slightly different vertical wavenumber than its neighbors when \(k_y = 0\) [25], [42]. As shown in Section 3 (see Figure 3), large-scale azimuthal motions can project onto a continuous spectrum of singular Alfvén eigenmodes and produce fine horizontal scales via phase mixing.

This effect is unimportant in the \(k_y = 0\) case [33] since the IGWs forced at the top of the domain are decoupled from azimuthal motions \(v\). However, in the \(k_y \neq 0\) case considered below, the IGWs have nonzero \(v\). To see this, first observe that the perturbations in density (\(\rho_0 \neq 0\)) associated with IGWs must be accompanied by pressure fluctuations (the \(z\) component of the momentum equation, 24 , implies that \(p_0 \neq 0\) if \(\rho_0 \neq 0\)). Then, the \(y\) components of the induction and momentum equations (combined as Equation 15 ) imply that these pressure fluctuations induce azimuthal motions (\(v_0 \neq 0\)) if \(k_y \neq 0\). The consequent coupling of the discrete modes to the Alfvén continuum results in the conversion to AWs and to a mixed SM-AW mode as discussed below.

4.2 Non-axisymmetric modes↩︎

The AWs and mixed SM-AW mode have features with fine scales in \(x\) associated with the singularity in (15 ), as shown later in this section. To regularize these features with minimal influence on the remaining large-scale modes, we introduce anisotropic diffusion to the induction equation: \[\label{eqn:anisoinduction} \partial_t \boldsymbol{b} = \boldsymbol{\nabla}\times(\boldsymbol{u}\times\boldsymbol{B}_0) + \eta_x \partial_x^2 \boldsymbol{b} + \eta (\partial_y^2 + \partial_z^2) \boldsymbol{b},\tag{18}\] with higher diffusivity in the \(x\) direction (\(\textit{Lu}_x = \omega l_z^2/\eta_x = 12.5\)) than in \(y\) or \(z\) (\(\textit{Lu} = 6.25 \times 10^4\)). We retain the leading order anisotropic diffusion terms in the reduced induction equation: \[\label{eqn:anisoinduction-reduced} -i \omega \boldsymbol{b}_{0 H} = i k_z B_{0 z} \boldsymbol{u}_{0 H} + \eta_x \partial_x^2 \boldsymbol{b}_{0 H} - \eta (k_y^2 + k_z^2) \boldsymbol{b}_{0 H},\tag{19}\] where \(\boldsymbol{u}_{0 H} = u_0 \boldsymbol{e}_x + v_0 \boldsymbol{e}_y\) and \(\boldsymbol{b}_{0 H} = b_{0 x}\boldsymbol{e}_x + b_{0 y} \boldsymbol{e}_y\). (Terms involving \(w_0\) and \(b_{0 z}\) are \(O(\textit{Fr})\) smaller and thus neglected). The full set of (resistive) reduced equations are derived in Appendix 7. Appendix 8.2 presents the reduced system as a generalized eigenvalue problem and details the procedure for finding the non-axisymmetric (\(k_y = 2\pi/L\)) resistive modes.

Figure 5b plots the real part of \(k_z\) for the \(k_y = 2\pi/L\) resistive modes. In the presence of diffusion, the ideal Alfvén continuum is replaced by discrete wavenumber branches [44]. However, these branches are so closely spaced for \(\textit{Lu}_x =12.5\) and \(\textit{Lu} = 6.25 \times 10^4\) that we plot the ideal Alfvén continuum (gray shaded region) for simplicity. As in the \(k_y = 0\) case, IGW-0 and IGW-1 (yellow and orange curves, respectively) closely follow the pure IGW wavenumber (gray vertical dash-dotted line) near the top of the domain. The two parities exhibit different behavior as they interact with background magnetic field.

There is a turning point at \(z/L \approx 0.0461\), where IGW-1 meets the SM-1 branch (purple curve in Figure 5b) and leaves behind an evanescent tail (evan.-1, dotted orange curve), as before. However, beyond the Alfvén wavenumber boundary (i.e., when \(\Re\{k_z (z)\} \ge k_{A B}(z)\), see Equation 16 ), SM-1 becomes a mixed SM-AW mode (SM-AW-1, dashed purple curve). The \(v_0\) eigenmodes for SM-AW-1 resemble SM-1 apart from sharp features at critical latitudes \(x_c\) where the mode resonates with the Alfvén continuum. Despite these resonances, SM-AW-1 behaves similarly to SM-1, refracting as it approaches the \(k_y = 0\) SM cutoff height (horizontal gray dashed line). The increase of the imaginary part of \(k_z\) with height (see Figure 5c, dashed purple curve) shows that SM-AW-1 decays faster than exponentially as it approaches the cutoff height.

This behavior is consistent with the results of [34] in spherical geometry. For non-axisymmetric dipole modes (\(\ell = \lvert m \rvert = 1\)), they identify a turning point between distinct IGW and SM branches. Further, they find that the SM eigenmodes develop sharp features at critical latitudes when the associated radial wavenumber satisfies the Alfvén dispersion relation somewhere in the domain. In this work, we characterize these modes as mixed SM-AW modes and distinguish them from AWs that do not directly connect to SM branches.

The mixed SM-AWs have large-scale features away from the critical latitudes and may be resolved with vertical diffusion (\(\eta \partial_z^2 \boldsymbol{b} \rightarrow -\eta k_z^2 \boldsymbol{b}_0\)) alone, as found by [30] and [34]. In contrast, the AWs are dominated by fine scale oscillations and are only resolved when the horizontal diffusion operator (\(\eta_x \partial_x^2 \boldsymbol{b}_0\)) is explicitly included in the WKB induction equation (26 ).

These AWs have wavenumbers with large imaginary part for the values of the diffusivities used here. Note that \(\Im\{k_z\}\) depends mainly on \(\eta_x\) here, since \(\eta_x \gg \eta\); we expect \(\Im\{k_z\}/\Re\{k_z\}\) would be much smaller if \(\textit{Lu}_x = \textit{Lu} = 6.25 \times 10^4\). Near the cutoff height, SM-AW-1 connects to one such branch of evanescent AWs (AW-1, dotted purple curve). Unlike SM-AW-1, the azimuthal velocity eigenmode \(v_0\) for AW-1 grows faster with height than \(u_0\); the AW transfers energy from latitudinal motions to azimuthal ones (though these still ultimately damp since \(v\) depends on the product \(v_0 \exp(i k_z z)\) where \(\Im\{k_z\} > 0\)).

The wave interactions for the opposite parity are slightly different. At \(z/L \approx 0.0522\), IGW-0 refracts into a sine-parity SM branch (SM-0, green curve in Figure 5b) and leaves behind an evanescent tail (evan.-0, dotted yellow curve). SM-0 refracts further as it propagates upwards, following the Alfvén wavenumber boundary \(k_z \approx k_{A B}(z)\). Near \(z/L \approx 0.163\), the SM-0 branch enters the Alfvén continuum (i.e., \(\Re\{k_z(z)\} > k_{A B}(z)\)) and connects to an evanescent AW branch (AW-0, dotted green curve).

4.3 Comparison between WKB theory and simulation↩︎

The reduced eigenvalue problem suggests that IGWs of both parity convert to decaying magnetohydrodynamic waves. To confirm this, we compare a WKB solution composed of the nine dipolar modes discussed above (IGW-0, IGW-1, evan.-0, evan.-1, SM-0, SM-1, SM-AW-1, AW-0, AW-1) to a third numerical simulation, IVP III, that uses anisotropic magnetic diffusion. In IVP III, we solve (18 ) instead of (5 ), with \(\textit{Lu}_x = 12.5\) and \(\textit{Lu} = 6.25 \times 10^4\) to match the non-axisymmetric eigenvalue problem.

Figure 6 compares the latitudinal velocity component of IVP III in the equilibrated state (\(u_{\text{IVP}}\), panel a) to the WKB solution (\(u_{\text{WKB}}\), panel b). A supplementary video temporarily available at this link shows the transient behavior of IVP III (details under Supplementary Materials). Because the WKB solution is only determined up to an arbitrary phase and amplitude, comparing it to the IVP solution requires a few steps, the full details of which are in Appendix 9. The overall amplitude and phase of \(u_{\text{WKB}}\) is determined by normalizing the WKB solution to \(u_{\text{IVP}}\) at an arbitrary point well above the cutoff heights and away from the zeros of the eigenmodes: \((x/L, z/L) = (1/3, 0.175)\). Crucially, we assume that there is total conversion from IGWs to SM waves at the turning points \(z_{t,\text{IGW-0}} \approx 0.0522 L\) and \(z_{t,\text{IGW-1}} \approx 0.0461 L\) for the sine- and cosine-parity waves, respectively. Further, we assume that the difference in phase \(\theta\) across each turning point is \(-\pi/2\), as in the axisymmetric case [33]. This assumption, along with that of total conversion between IGWs and SM waves, is supported by the excellent match between \(u_{\text{IVP}}\) and \(u_{\text{WKB}}\) outside of the damping regions (hatched areas) and below the forcing height (black dotted line).

To enable detailed comparison between the simulation and WKB theory, Figure 6c,d shows the sine-parity components of \(u_{\text{IVP}}\) and \(u_{\text{WKB}}\), the latter consisting only of IGW-0, evan.-0, SM-0, and AW-0 modes. The numerical simulation agrees with the results of the WKB analysis assuming total conversion from IGW-0 to SM-0 at the turning point \(z = z_{t,\text{IGW-0}}\) (yellow line). SM-0 then damps as it approaches a magnetic cutoff height near \(z = z_{t,\text{AW-0}} \approx 0.163 L\) (the diamond-shaped features in Figure 6c,d vanish near the dotted green line).

Figure 6e,f shows just the cosine-parity components of \(u_{\text{IVP}}\) and \(u_{\text{WKB}}\), the latter consisting only of IGW-1, evan.-1, SM-1, SM-AW-1, and AW-1 modes. Again, the match between numerical simulation and WKB theory supports total conversion from IGW to SM wave at the turning point \(z=z_{t,\text{IGW-1}}\) (orange line). SM-1 then crosses the Alfvén wavenumber boundary (at the dashed purple line) and completely converts to a mixed SM-AW mode (SM-AW-1). SM-AW-1 damps rapidly as it approaches \(z = z_{t,\text{AW-1}} \approx 0.091 L\) (the lenticular features in Figure 6e,f vanish near the dotted purple line).

Figure 7 compares the azimuthal velocity component of IVP III in the equilibrated state (\(v_{\text{IVP}}\), panel a) to the WKB solution (\(v_{\text{WKB}}\), panel b). The latter is constructed using the same amplitude and phase shifts as \(u_{\text{WKB}}\). The agreement between WKB theory and simulation is not as good for \(v\) as it is for \(u\), owing to the strong influence of the Alfvén continuum on \(v\) (discussed below).

Figure 7c,d shows the cosine components of \(v_{\text{IVP}}\) and \(v_{\text{WKB}}\). Above the turning point (yellow line), SM-0 develops ridge-like features centered at \(x/L = 0,0.5,1\) that sharpen as it approaches the Alfvén continuum in wavenumber space. These ridges lie at the critical latitudes associated with the Alfvén boundary \(k_z(z) = k_{A B}(z)\) (vertical black dashed lines).

Near \(z = z_{t,\text{AW-0}}\) (Figure 7c,d, dotted green line), SM-0 crosses into the Alfvén continuum and converts to an evanescent AW (AW-0, faint chevron-shaped features above the dotted green line). Comparison to IVP I (Figure 2e) and IVP II (Figure 4e) shows that the damping rate of the large chevron-shaped features associated with AW-0 depends strongly on the diffusivities. Though phase-mixing is present in IVP I and II, it acts less efficiently on these features since they occur where the local Alfvén wavenumber varies most slowly (\(x/L = 0, 0.5, 1\)). Regardless, the horizontal structure of these AWs is distinct from that of the IGW eigenmodes, and they are not expected to couple to up-going IGWs even in the limit of vanishing diffusivity.

Figure 7e,f plots the sine components of \(v_{\text{IVP}}\) and \(v_{\text{WKB}}\). SM-AW-1 develops sharp shingle-like features at critical latitudes (black curves) where the mode resonates with the Alfvén continuum (i.e., where \(k_z^2 = \omega^2/v_{A z}^2\)). Without resistivity, these sharp features correspond to singularities in the azimuthal velocity eigenmode \(v_0\). Anisotropic diffusion regularizes these peaks, though it acts slightly differently in the eigenproblem than in the simulation. As a result, the shingles are sharper and have higher amplitude in \(v_{\text{WKB}}\) than in \(v_{\text{IVP}}\).

It is not entirely surprising that we miss some details of the diffusion-mediated behavior, as it is difficult to include diffusion terms self-consistently in an asymptotic analysis [45]. While we have tried different ways of including diffusion into our asymptotic model, we have not been able to recover the precise behavior of these diffusive modes seen in the IVP. We were only able to obtain eigenfunctions that share the the IVP’s broad, smooth peaks at SM-AW-1’s critical latitudes by increasing \(\eta_x\) such that \(\textit{Lu}_x=12.5\). This value of \(\eta_x\) is significantly larger than the horizontal diffusivity needed to produce similar features in IVP I (\(\textit{Lu}_x = \textit{Lu}= 6.25 \times 10^4\)).

Regardless of these differences, the shingle-like features in both IVP III and the WKB solution damp after SM-AW-1 converts to AW-1 (Figure 7e,f, above the purple dotted line). For lower horizontal diffusivity, energy in SM-AW-1 scatters among AWs that undergo rapid phase mixing with height, as seen in IVP I (Figure 2f) and IVP II (Figure 4f). These AWs damp, and no up-going large-scale features reach the top of the domain.

The behavior of the large-scale wave modes discussed above is supported by the match between IVP III and the WKB solution. To better compare the simulation and theory, Figure 8 overplots vertical profiles of IVP III (blue dashed curve) and the WKB solution (red curve) at several different latitudes \(x/L = 1/3,0.5, 0.75, 0.95\). Panels a and e plot \(u\) and \(v\), respectively, at the fitting latitude \(x/L = 1/3\). The only major differences between theory and simulation lie above the forcing height (black dotted line) and at the turning points (yellow and orange lines), where the WKB solution is expected to diverge. The excellent agreement between \(u_{\text{IVP}}\) and \(u_{\text{WKB}}\) persists at \(x/L = 0.5, 0.75, 0.95\) (Figure 8b,c,d).

As discussed above, the quality of fit is slightly diminished for the azimuthal flow \(v\). Figure 8f plots \(v\) at \(x/L = 0.5\), i.e., only the modes for which \(v_0\) has cosine parity (IGW-0, evan.-0, SM-0, AW-0). The ridge-like features of SM-0 (see Figure 7c,d) are sharper in the WKB theory than in the IVP; \(v_{\text{WKB}}\) has higher amplitude than \(v_{\text{IVP}}\) between the yellow and dotted green lines in panel f. The WKB solution also overestimates \(v\) at \(x/L = 0.95\) (panel h) between the Alfvén boundary intersection height (dashed purple line) and the SM-AW-1 cutoff height (dotted purple line). This mismatch is due to the shingle-like features in SM-AW-1 discussed above (see Figure 7e,f) that are sharper in the theory than in the simulation.

These differences notwithstanding, the overall agreement between the WKB solution and simulation show that down-going IGWs forced at the top of the domain convert completely to a mix of SM waves and AWs, which develop fine scales and damp. Although we could only explicitly compare the WKB theory to the simulation with anisotropic diffusion (IVP III), the same behavior is observed in the large-scale modes in IVP I and IVP II. Thus, all three simulations support this interpretation.

5 Conclusions↩︎

This work analyzes the interaction of non-axisymmetric IGWs with spatially-varying magnetic fields using numerical simulations. We extend the Cartesian model of [33] into the third dimension (\(y\), azimuth) such that the IGWs and SM waves in their model are accompanied by a continuous spectrum of singular AW modes, as found in the WKB analyses of [30] and [34] in spherical geometry. While these works are geometrically accurate, they investigate only the eigenvalues and eigenmodes of the reduced equations. By staying in Cartesian geometry, we are able to go beyond [34] and obtain an amplitude equation, allowing us to compare a WKB solution directly to numerical simulations and to study the conversion of modes across turning points.

The agreement between our simulations and the WKB theory shows that down-going non-axisymmetric IGWs convert completely to up-going SM waves (as in the axisymmetric case, [33]). The SM waves then interact with the Alfvén continuum in different ways, depending on their parity. The SM wave for which the azimuthal flow vanishes where the vertical background magnetic field component is strongest (SM-1) is at first weakly influenced by the Alfvén continuum. SM-1 converts to a mixed SM-AW mode (SM-AW-1), which resembles SM-1 apart from critical layers located at latitudes where the real part of the vertical wavenumber satisfies the Alfvén dispersion relation. The SM wave of the opposite parity (SM-0, for which the peaks in azimuthal flow and the vertical component of the background field coincide) is more strongly influenced by the continuous spectrum. Instead of refracting into the Alfvén continuum like SM-1, SM-0 propagates higher with a vertical wavenumber that follows the Alfvén wavenumber boundary \(k_{A B}(z)\) (16 ). Regardless, both SM-AW-1 and SM-0 approach magnetic cutoff heights near which the latitudinal flow rapidly damps. These results agree with the predictions of [30] and [34], who identified similar modes in spherical eigenvalue problems.

At the magnetic cutoff heights, the SM-AW and SM wave convert to AWs. Phase mixing of AWs produces fine scales that result in efficient damping in the presence of diffusion. Simulations with different magnetic diffusivities (\(\textit{Lu} = 6.25 \times 10^4\) in IVP I; \(\textit{Lu} = 6.25 \times 10^5\) in IVP II) show that the lower the diffusivity, the further the AWs propagate. Despite these differences, the large-scale modes are the same and agree with the WKB solution in all simulations.

Our results have important implications for the asteroseismic signature of strongly-magnetized stars. In many massive stars, surface convection excites nonradial acoustic waves [46] that couple to IGWs in the radiative interior [47], [48]. Typically, these IGWs reflect off of the inner region of the core and set up global mixed modes [49] with predictable fluctuations in surface brightness [11]. However, in roughly 20% of RGB stars, dipole pulsation modes have lower-than-expected amplitudes [22], [23]. [33] explain the suppression of axisymmetric dipole modes via interactions of down-going IGWs with core magnetic fields above a critical strength.

Importantly, suppressed amplitudes are also observed for non-axisymmetric dipole modes, and there is no clear pattern of non-axisymmetric modes being suppressed more or less than axisymmetric modes. The present work shows that these observations could be explained by similar wave interactions with strong magnetic fields. Non-axisymmetric IGWs convert to a mix of SM and AWs that damp due to the production of fine scales via refraction (in the case of SM waves) and phase-mixing (in the case of AWs).

It has been suggested that these AWs could travel along closed field loops in the core and convert to up-going IGWs [34]. If so, an alternate mechanism would be required to explain the observed suppression of non-axisymmetric dipole modes. However, we find no evidence of conversion from AWs back to IGWs in our simulations. Though the diffusivities are far lower in RGB stars than in our simulations, we expect efficient damping of all AWs near the top of the core where the magnetic field should experience faster-than-exponential decay with radius. Thus, both axisymmetric and non-axisymmetric dipole mixed modes should have suppressed amplitudes in stars with sufficiently strong core magnetic fields, and our simulations do not provide a mechanism explaining the g-mode character of partially-suppressed mixed modes described in [29].

These results bolster recent efforts [30]–[32] to infer the strength of magnetic fields in astrophysical bodies based on wave suppression. We find that non-axisymmetric IGWs convert to SM waves where the magnetic field strength reaches only a fraction of the value needed for the conversion of axisymmetric IGWs (we observe that for the same background field strength, the IGW\(\rightarrow\)SM turning points are higher in the non-axisymmetric case; see Figure 5a,b). Thus, we expect the critical field strength deduced from the axisymmetric problem [33] to provide a good lower bound on the magnetic field needed to account for the suppression of both axisymmetric and non-axisymmetric dipole modes.

Future work should refine our results with numerical simulations in spherical geometry. The Cartesian model used here simplifies the problem and facilitates the derivation of the WKB amplitude equation, which is essential for comparing theory to simulation. Though the Cartesian wave modes discussed in this work largely resemble the spherical modes found by [34], some qualitative differences exist. For example, [34] find (what we interpret here as) mixed SM-AW modes with large-scale oscillations confined to an equatorial band between Alfvén critical latitudes. It is unclear how much this equatorial focusing effect changes the wave interaction problem, since the behavior of the SM-AW modes in Cartesian and spherical geometries is otherwise similar. Further investigations should also account for rotation, which acts differently on prograde versus retrograde modes (i.e., modes with \(m\) of opposite sign, [50]). Though rotational and geometric effects may introduce the differences described above, we expect our main conclusions to hold: energy carried inward by both axisymmetric and non-axisymmetric IGWs is ultimately lost via interactions with a strong magnetic field.

Supplementary Materials↩︎

Animated versions of Figures 2, 3, and 4 are temporarily available here. Additionally, a supplementary video showing the transient behavior of IVP III is temporarily available here; animations of \(u_{\text{IVP}}\) and \(v_{\text{IVP}}\) (decomposed by parity) are overlaid with horizontal lines (initially at the forcing height) that move vertically according to the group velocity of each WKB mode (computed using equation 30 in Appendix 7). The final frame corresponds to the equilibrated state in panels c and e of Figures 6 and 7.

6 Dimensionless equations↩︎

The numerical simulations and WKB analysis use a dimensionless form of the magneto-Boussinesq equations (1 ) that are derived as follows. Quantities denoted with tildes (\(\tilde{\;\;}\)) are dimensionless. We scale time as \(t = \omega^{-1} \tilde{t}\) and the background magnetic field as \(\boldsymbol{B}_0 = \mathcal{B} \tilde{\boldsymbol{v}}_A =\mathcal{B}(\tilde{v}_{A x} \boldsymbol{e}_x + \tilde{v}_{A z} \boldsymbol{e}_z)\). Based on the dispersion relation for pure IGWs (7 ), we anticipate large horizontal length scales (comparable to forcing wavelength \(L\)) and small vertical scales (comparable to \(l_z = L \omega /N\)). Accordingly, we scale distances anisotropically: \((x,y) = (L\tilde{x},L\tilde{y})\), \(z = l_z \tilde{z}\), and \(k_y = \tilde{k}_y/L\). The velocity and magnetic field perturbations are also scaled anisotropically to balance their respective divergence-free constraints (3 ,6 ). We decompose each as \(\boldsymbol{u} = \boldsymbol{u}_H + w\boldsymbol{e}_z\) and \(\boldsymbol{b} = \boldsymbol{b}_H + b_z\boldsymbol{e}_z\), and let \(\boldsymbol{u}_H = \mathcal{U} \tilde{\boldsymbol{u}}_H\), \(w = \mathcal{U} (l_z/L) \tilde{w}\), \(\boldsymbol{b}_H = \mathcal{U} \sqrt{\mu_0 \overline{\rho}} \tilde{\boldsymbol{b}}_H\), and \(b_z = \mathcal{U} (l_z/L) \sqrt{\mu_0 \overline{\rho}} \tilde{b}_z\), where \(\mathcal{U}\) is the (velocity) amplitude of the IGWs. The mass equation (4 ) is balanced by scaling the density perturbation as \(\rho = \overline{\rho} N^2 l_z \mathcal{U}/(g \omega L)\tilde{\rho}=(\bar \rho N {\cal U}/g) \tilde{\rho}\). Finally, the pressure perturbation is scaled to reflect leading-order hydrostatic balance in the vertical momentum equation: \(p = \overline{\rho}N^2 l_z^2\mathcal{U}/(\omega L) \tilde{p}= (\bar \rho \omega\mathcal{U}L)\tilde{p}\). After nondimensionalizing and dropping tildes (\(\tilde{\;\;}\)), the linearized equations become \[\label{eqn:ndimmagbouss} \begin{equation} \partial_t \boldsymbol{u}_H + \boldsymbol{\nabla}_H p - \Gamma v_{A z} \partial_z \boldsymbol{b}_H = -\textit{Fr}^2\Gamma v_{A z} \boldsymbol{\nabla}_H b_z + \textit{Fr} \Gamma (\boldsymbol{\nabla}_H \times \boldsymbol{b}_H) \times v_{A x} \boldsymbol{e}_x, \end{equation} \begin{equation} \partial_z p + \rho= -\textit{Fr}^2\partial_t w + \textit{Fr} \Gamma v_{A x}(\textit{Fr}^2\partial_x b_z - \partial_z b_x), \end{equation} \begin{equation} \boldsymbol{\nabla}_H \boldsymbol{\cdot} \boldsymbol{u}_H + \partial_z w = 0, \end{equation} \begin{equation} \partial_t \rho - w = 0, \end{equation} \begin{multline} \partial_t \boldsymbol{b}_H - \Gamma v_{A z} \partial_z \boldsymbol{u}_H - \textit{Lu}_x^{-1} \textit{Fr}^2\partial_x^2\boldsymbol{b}_H - \textit{Lu}^{-1}\left(\textit{Fr}^2\partial_y^2 +\partial_z^2\right)\boldsymbol{b}_H = -\textit{Fr}^2\Gamma w \partial_\mathcal{Z}v_{A x}\boldsymbol{e}_x\\ +\textit{Fr} \Gamma \left(\boldsymbol{v}_A \boldsymbol{\cdot} \boldsymbol{\nabla}_H \boldsymbol{u}_H -\boldsymbol{u}_H \boldsymbol{\cdot} \boldsymbol{\nabla}_H v_{A x}\boldsymbol{e}_x\right), \end{multline} \begin{equation} \partial_t b_z - \Gamma(v_{A z}\partial_z w - u\partial_x v_{A z}) - \textit{Lu}_x^{-1} \textit{Fr}^2\partial_x^2 b_z - \textit{Lu}^{-1}\left(\textit{Fr}^2\partial_y^2 + \partial_z^2\right) b_z= \textit{Fr} \Gamma (v_{A x}\partial_x w - w \partial_\mathcal{Z}v_{A z}), \end{equation} \begin{equation} \boldsymbol{\nabla}_H \boldsymbol{\cdot} \boldsymbol{b}_H + \partial_z b_z = 0, \end{equation}\tag{20}\] where \(\boldsymbol{\nabla}_H = \boldsymbol{e}_x \partial_x + \boldsymbol{e}_y \partial_y\), and where we have assumed that the background field is current-free (i.e., \(\partial_{\mathcal{Z}} v_{A x} - \partial_x v_{A z} = 0\)). For anisotropic diffusion (IVP III), \(\textit{Lu}_x \neq \textit{Lu}\). For the remainder of the Appendix, all quantities are dimensionless and will be written without tildes (\(\tilde{\;\;}\)).

Finally, we observe that, in the absence of diffusion, the horizontally-integrated wave energy is governed by \[\begin{gather} \label{eqn:horizavgenergy} {\partial_t}\int_{-\pi/k_y}^{-\pi/k_y} \int_0^1 \frac{1}{2}\left(\boldsymbol{u}_H \cdot \boldsymbol{u}_H + \textit{Fr}^2 w^2 + {\rho^2} + \boldsymbol{b}_H \cdot \boldsymbol{b}_H + \textit{Fr}^2 b_z^2\right) \mathrm{d}x\mathrm{d}y \\ + \partial_z \int_{-\pi/k_y}^{-\pi/k_y} \int_0^1 \left[p w - \Gamma(\boldsymbol{b}_H \cdot \boldsymbol{u}_H) v_{A z} + \Gamma \textit{Fr} b_x v_{A x} w \right] \mathrm{d}x\mathrm{d}y= 0, \end{gather}\tag{21}\] which is used to derive an expression for the WKB group velocity in the following section.

7 Development of the reduced equations↩︎

We set \(\epsilon = \textit{Fr}\) and expand each field \(u\), \(v\), \(w\), \(p\), \(\rho\), \(b_x\), \(b_y\), \(b_z\) as, e.g., \[p = \Re\left\{\hat{p}(x, \mathcal{Z})\exp\left(i \int k_z(\mathcal{Z})\mathrm{d}z +i k_y y - i t\right)\right\},\] with \[\label{eqn:wkbexpansion} \hat{p}(x,\mathcal{Z}) = A(\mathcal{Z}) p_0(x, \mathcal{Z}) + \epsilon p_1(x, \mathcal{Z}) + \epsilon^2 p_2(x, \mathcal{Z}) +...\tag{22}\] where the amplitude function \(A\), vertical wavenumber \(k_z\), and each \(p_i\) vary with the “slow” vertical coordinate \(\mathcal{Z} = \epsilon z\). Note that \(k_z\) here is nondimensionalized by \(l_z\). Then, the derivative operators transform as, e.g., \(\partial_t p = -i \hat{p} \exp(i \Phi)\), \(\partial_x p = \partial_x \hat{p} \exp(i \Phi)\), \(\partial_y p = i k_y \hat{p} \exp(i \Phi)\), and finally, \(\partial_z p = i k_z({\cal Z}) \hat{p} \exp(i \Phi) + \epsilon \partial_{\mathcal{Z}} \hat{p} \exp(i \Phi)\), where \(\Phi = \int k_z(\mathcal{Z})\mathrm{d}z + k_y y - t\).

We substitute these ansätze into (20 ) and analyze the distinguished limit \(\Gamma = O(1)\) as \(\epsilon \to 0\). From the left-hand side of (20 ) we deduce the leading order resistive equations, \[\tag{23} \begin{equation} -i \boldsymbol{u}_{0 H} + \hat{\boldsymbol{\nabla}} p_0 - i \Gamma v_{A z}k_z\boldsymbol{b}_{0 H} = \boldsymbol{0}, \end{equation} \begin{equation}\tag{24} i k_z p_0 + \rho_0 = 0, \end{equation} \begin{equation} \hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} \boldsymbol{u}_{0 H} + i k_z w_0 = 0, \end{equation} \begin{equation}\tag{25} -i \rho_0 - w_0= 0, \end{equation} \begin{equation}\tag{26} -i \boldsymbol{b}_{0 H} - i \Gamma k_z v_{A z} \boldsymbol{u}_{0 H} -\textit{Lu}_x^{-1}\textit{Fr}^2\partial_x^2 \boldsymbol{b}_{0 H} + \textit{Lu}^{-1}(\textit{Fr}^2 k_y^2+k_z^2)\boldsymbol{b}_{0 H}= \boldsymbol{0}, \end{equation}\] where \(\hat{\boldsymbol{\nabla}} = \boldsymbol{e}_x \partial_x + \boldsymbol{e}_y i k_y\), \(\boldsymbol{u}_{0 H} = u_0 \boldsymbol{e}_x + v_0 \boldsymbol{e}_y\), and \(\boldsymbol{b}_{0 H} = b_{0 x} \boldsymbol{e}_x + b_{0 y} \boldsymbol{e}_y\). Note that, while the diffusion terms in (23 ) are formally higher order in \(\textit{Fr}\), we retain them in the equations for regularization.

In the absence of diffusion (\(\textit{Lu}_x^{-1} = \textit{Lu}^{-1} = 0\)), the leading order system (23 ) can be combined into a single homogeneous equation, \(L_{k_z}p_0 = 0\) where \[L_{k_z}(\cdot) = k_z^2 (\cdot) + \hat{\boldsymbol{\nabla}} \boldsymbol{\boldsymbol{\cdot}} \left( \frac{\hat{\boldsymbol{\nabla}} (\cdot)}{1-\Gamma^2 v_{A z}^2 k_z^2 } \right).\]

The remaining piece of the leading-order WKB solution is the amplitude function \(A(\mathcal{Z})\), which is obtained at the next order. The \(O(\epsilon)\) equations form an inhomogeneous linear system that may be combined into the single equation \[\label{eqn:ordereps} L_{k_z}p_1 = f,\tag{27}\] where \[\begin{equation} f = -i\partial_{\mathcal{Z}}(A w_ 0) +i {k_z}\partial_{\mathcal{Z}} (A p_ 0) - {k_z^2}\Gamma A v_{A x} b_{0 x} + \hat{\boldsymbol{\nabla}} \cdot \left( \frac{\boldsymbol{h}}{1-\Gamma^2 k_z^2 v_ {A z}^2}\right), \end{equation} \begin{equation} \boldsymbol{h} =\Gamma v_ {A z} \partial_{\mathcal{Z}}(A \boldsymbol{b}_{0 H} )+ \Gamma ( \hat{\boldsymbol{\nabla}} \times A \boldsymbol{b}_{0 H})\times v_ {A x}\boldsymbol{e}_x - {\Gamma^2 v_ {A z} k_z} \left[v_ {A z} \partial_{\mathcal{Z}} (A \boldsymbol{u}_{0 H}) + v_ {A x} \partial_x(A \boldsymbol{u}_{0 H}) - A u_ 0 \partial_x v_ {A x} \boldsymbol{e}_x \right]. \end{equation}\] A solvability condition on the inhomogeneous term \(f\) is required to ensure that a bounded solution exists, and may be obtained by taking the inner product of (27 ) with an element of the kernel of the adjoint operator \(L_{k_z}^\dagger\). We observe that \(L_{k_z}\) is self-adjoint up to complex conjugation of \(k_z\) (i.e., \(L_{k_z}^\dagger = L_{k_z^\ast}\)) for the inner product \(\langle a,b\rangle = \int_0^1 a^\ast b \mathrm{d}x\). Using this fact, we have that \(p_0^\ast\) is in the kernel of \(L_{k_z}^\dagger\); that is, \(L_{k_z}^\dagger p_0^\ast = L_{k_z^\ast}p_0^\ast =(L_{k_z}p_0)^\ast = 0\). Thus, taking the inner product of (27 ) with \(p_0^\ast\) yields the desired solvability condition: \[\label{eqn:solvcond} 0 = \langle p_0^\ast, f \rangle .\tag{28}\] Equation (28 ) simplifies to a first order differential equation in \(\mathcal{Z}\) which determines the amplitude function up to a multiplicative constant: \[\label{eqn:amp} A(\mathcal{Z}) = \left\{\int_0^1 \left[p_0 w_0 + \Gamma v_{A z} (b_{0 x} u_0 - b_{0 y} v_0)\right]\mathrm{d}x \right\}^{-1/2}.\tag{29}\] A detailed derivation of (29 ) may be found in [51].

Finally, if the amplitude \(A\) is allowed to vary slowly with time \(\mathcal{T} = \epsilon t\) (before the solution equilibrates), then one may show (by substituting the WKB ansätze into Equation 21 and keeping only the leading order terms) that the WKB wave energy \[E_0 =\frac{\lvert A(\mathcal{Z},\mathcal{T}) \rvert^2}{4}\int_0^1 \left( {\boldsymbol{u}_{0 H}\cdot\boldsymbol{u}_{0 H}^\ast} + {\rho_0 \rho_0^\ast} + {\boldsymbol{b}_{0 H}\cdot\boldsymbol{b}_{0 H}^\ast}\right) \mathrm{d}x\] is transported vertically (in the absence of diffusion) according to \[(\partial_{\mathcal{T}} + c_{g z} \partial_\mathcal{Z}) E_0 = -E_0 \partial_\mathcal{Z} c_{g z},\] where the vertical WKB group velocity is given by \[\label{eqn:groupvel} c_{g z} = \frac{\text{(rate of) pressure work} + \text{Poynting flux}}{\text{wave energy}} = \frac{\int_0^1 \Re\left\{p_0^\ast {w}_ 0- \Gamma v_{A z} {\boldsymbol{b}}_ {0 H} \cdot {\; \boldsymbol{u}}_{0 H}^\ast\right\} \mathrm{d}x}{\frac{1}{2}\int_0^1 \left( {\boldsymbol{u}_{0 H}\cdot\boldsymbol{u}_{0 H}^\ast} + {\rho_0 \rho_0^\ast} + {\boldsymbol{b}_{0 H}\cdot\boldsymbol{b}_{0 H}^\ast}\right) \mathrm{d}x}.\tag{30}\]

8 Eigenvalue problem implementation in Dedalus↩︎

8.1 Ideal axisymmetric eigenproblem↩︎

The WKB wave modes and vertical wavenumbers are found numerically using the Dedalus pseudospectral code. For the ideal “axisymmetric” case (\(\textit{Lu}_x^{-1} = \textit{Lu}^{-1} = 0\), \(k_y =0\)) in Figure 5a, we solve the following generalized eigenvalue problem \[\label{eqn:idealevp} \begin{bmatrix} 0 & i\partial_x & - k_y \\ i \partial_x & 1 & 0 & \\ - k_y & 0 & 1 & \end{bmatrix} \begin{bmatrix} p_0 \\ u_0 \\ v_0 \end{bmatrix} +k_z^2 \begin{bmatrix} 1 & 0 & 0 \\ 0&-\Gamma^2 v_{A z}^2 & 0 \\ 0&0 & -\Gamma^2 v_{A z}^2 \end{bmatrix} \begin{bmatrix} p_0 \\ u_0 \\ v_0 \end{bmatrix} = \boldsymbol{0}\tag{31}\] for eigenvalue \(k_z^2\) and eigenvector \([p_0 \;u_0 \; v_0 ]^T\). The system (31 ) only depends on \(\mathcal{Z}\) parametrically (through \(v_{A z}\)) and may be solved over one dimension (\(x\)) at each level in \(\mathcal{Z}\). We use Dedalus’ dense solver to find the eigenpairs using \(N_x = 96\) complex Fourier modes in \(x\) at each of \(N_z = 512\) uniformly spaced values of \(\mathcal{Z} \in [0,0.25]\). To assess numerical convergence, the process is repeated with \(N_x = 128\) complex Fourier modes, and only the eigenvalues that differ by less than \(10^{-8}\) times their original value are kept.

8.2 Resistive non-axisymmetric eigenproblem↩︎

A different procedure yields the resistive wave modes used to construct the non-axisymmetric (\(k_y=2\pi\)) WKB solution in Section 4. We find the generalized eigenvalues \(k_z\) (see Figure 5b,c) and eigenvectors \(\boldsymbol{q}_0 = [u_0 \;v_0\; w_0\;p_0\;\rho_0\;b_{0 x}\;b_{0 y}\;b_{0 x,z}\;b_{0 y,z}]^T\) of \[\label{eqn:resistiveevp} \begin{bmatrix} 1 & 0 & 0 & i\partial_x & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & -k_y & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0\\ i\partial_x & -k_y & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 &1 -i\hat{D} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 -i\hat{D} & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ \end{bmatrix} \begin{bmatrix} u_0 \\v_0\\ w_0\\p_0\\b_{0 x}\\b_{0 y}\\b_{0 x,z}\\b_{0 y,z} \end{bmatrix} + k_z \begin{bmatrix} 0 & 0 & 0 & 0 & \Gamma v_{A z} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \Gamma v_{A z} & 0 & 0\\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0\\ \Gamma v_{A z} & 0 & 0 & 0 & 0 & 0 & \textit{Lu}^{-1} & 0\\ 0 & \Gamma v_{A z} & 0 & 0 & 0 & 0 & 0 & \textit{Lu}^{-1}\\ 0 & 0 & 0 & 0 & -i & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & -i & 0 & 0\\ \end{bmatrix} \begin{bmatrix} u_0 \\v_0\\ w_0\\p_0\\b_{0 x}\\b_{0 y}\\b_{0 x,z}\\b_{0 y,z} \end{bmatrix} = \boldsymbol{0}\tag{32}\] where \(\hat{D} = \textit{Lu}_x^{-1}\textit{Fr}^2\partial_x^2 - \textit{Lu}^{-1}\textit{Fr}^2k_y^2\) is the transformed horizontal diffusion operator. The system (32 ) is formed by eliminating \(\rho_0\) from (23 ) and appending the equations \(b_{0 x,z} = i k_z b_{0 x}\) and \(b_{0 y,z} = i k_z b_{0 y}\).

In the presence of diffusion, (32 ) yields four branches of “dipolar” wavenumbers \(k_z\) (i.e., those for which \(\lvert p_0\rvert\) has two zeros on \(x \in [0,1)\)) that vary continuously over \(\mathcal{Z} \in [0,0.25]\). Our procedure for finding the eigenpairs along these branches is as follows. First, we set \(\mathcal{Z}=0\) in (32 ) and apply Dedalus’ dense solver at two different resolutions (\(N_x = 96\) and \(N_x=128\) complex Fourier modes in \(x\)). At this height, all four eigenvalues correspond to evanescent waves. Next, the resolved eigenvalue for each of the four dipolar branches is used as the target for Dedalus’ sparse solver with \(N_x = 1024\) complex Fourier modes. Then, we increment \(\mathcal{Z}\) by \(\Delta \mathcal{Z} = 0.25/N_z\) (with \(N_z = 3072\)) and use the eigenvalue from the previous level in \(\mathcal{Z}\) as the target for the sparse solver at the current level. The process is repeated up to \(\mathcal{Z} = 0.25\), keeping only the mode with the same parity as the previous level. Thus, at every \(j\)th level, we keep one eigenpair with eigenvalue \(k_{z,j}\) and eigenfunctions \(\boldsymbol{q}_{0,j} = [u_{0,j}(x) \;\;v_{0,j}(x)\;\; w_{0,j}(x)\;\;p_{0,j}(x) \dots]^T\). At this stage, the set of eigenfunctions \(\boldsymbol{q}_{0,j}\) for each of the four branches is normalized such that \(p_{0,j}(x=1/3)\) is constant for all \(j=1,2,\dots, N_z\) (this ensures continuity over \(\mathcal{Z}\)).

The four wavenumber branches have turning points \(\mathcal{Z}_t\) (which we identify as maxima in \(\lvert \mathrm{d} k_z/ \mathrm{d} \mathcal{Z}\rvert\)), where the assumption that \(\mathrm{d} k_z/ \mathrm{d} \mathcal{Z} = O(\epsilon)\) is violated. Thus, the WKB solution is only asymptotically exact away from the turning points, and these heights naturally delimit different WKB wave modes (e.g., a turning point separates the sine-parity IGW mode, IGW-0, from the associated evanescent mode, evan.-0). In addition, the cosine-parity SM mode (SM-1) is separated from the mixed SM-AW mode (SM-AW-1) at the height \(\mathcal{Z}_{A B}\) at which their vertical wavenumbers intersect the Alfvén wavenumber boundary, \(k_z = k_{A B}\) (see Figure 5b, where solid purple and dashed purple curves intersect the edge of the gray shaded region).

The turning points \(\mathcal{Z}_t\) and Alfvén boundary intersection height \(\mathcal{Z}_{A B}\) split the four continuous branches into 11 separate wave modes. Two of these modes have \(\Im\{k_z\} > 0\) at \(\mathcal{Z} =0\) such that the waves grow exponentially with depth near the bottom of the domain. We discard these two unphysical evanescent modes, leaving the nine wave modes that comprise the WKB solution in Section 4. The first branch is split into evan.-0 (for \(0 \le \mathcal{Z} < \mathcal{Z}_{t,\text{IGW-0}} \approx 0.0522\)) and IGW-0 (for \(\mathcal{Z} \ge \mathcal{Z}_{t,\text{IGW-0}}\)); the second branch is split into evan.-1 (for \(0 \le \mathcal{Z} < \mathcal{Z}_{t,\text{IGW-1}} \approx 0.0461\)) and IGW-1 (for \(\mathcal{Z} \ge \mathcal{Z}_{t,\text{IGW-1}}\)); the third branch is split into SM-0 (for \(\mathcal{Z}_{t,\text{IGW-0}} \le \mathcal{Z} < \mathcal{Z}_{t,\text{AW-0}} \approx 0.163\)) and AW-0 (for \(\mathcal{Z} \ge \mathcal{Z}_{t,\text{AW-0}}\)); and the fourth branch is split into SM-1 (for \(\mathcal{Z}_{t,\text{IGW-1}} \le \mathcal{Z} < \mathcal{Z}_{A B} \approx 0.0496\)), SM-AW-1 (for \(\mathcal{Z}_{A B} \le \mathcal{Z} < \mathcal{Z}_{t,\text{AW-1}} \approx 0.0910\)), and AW-1 (for \(\mathcal{Z} \ge \mathcal{Z}_{t,\text{AW-1}}\)).

For each of the nine waves, the full eigenmodes \(\boldsymbol{q}_0(x,\mathcal{Z}) = [u_0(x,\mathcal{Z}) \;\;v_0(x,\mathcal{Z})\;\; w_0(x,\mathcal{Z})\;\;p_0(x,\mathcal{Z}) \dots]^T\) are formed by “stitching” together the eigenfunctions \(\boldsymbol{q}_{0,j}(x) = [u_{0,j}(x) \;\;v_{0,j}(x)\;\; w_{0,j}(x)\;\;p_{0,j}(x) \dots]^T\) at each \(\mathcal{Z}_j\) over the vertical extent of the wave. The eigenmodes for IGW-0 and IGW-1 are normalized such that \(\langle p_0,p_0\rangle=1\) at the forcing height \(\mathcal{Z}_0\) with \(p_{0,\text{IGW-0}} = 1\) at \((x,\mathcal{Z})=(0,\mathcal{Z}_0)\) and \(p_{0,\text{IGW-1}} = 1\) at \((x,\mathcal{Z})=(1/4,\mathcal{Z}_0)\). This ensures that \(p_{0,\text{IGW-0}} + i p_{0,\text{IGW-1}} \approx \exp(i k_x x)\) at \(\mathcal{Z} = \mathcal{Z}_0\) to match the forcing. Then, the eigenmodes for SM-0 and SM-1 are rescaled such that \(p_{0,\text{SM-0}} \approx p_{0,\text{IGW-0}}\) at \((x,\mathcal{Z}) = (0,\mathcal{Z}_{t,\text{IGW-0}})\) and \(p_{0,\text{SM-1}} \approx p_{0,\text{IGW-1}}\) at \((x,\mathcal{Z)} = (1/4,\mathcal{Z}_{t,\text{IGW-1}})\). (Since we assume perfect conversion from IGWs to SM waves, the eigenmodes must match across the turning points). Finally, the eigenmodes for the remaining waves are rescaled such that for every pair of adjacent wave modes from the same branch (e.g., evan.-0 and IGW-0, but not evan.-0 and SM-0), \(p_0\) is continuous across the associated transition height (either a turning point or \(\mathcal{Z}_{A B}\)).

9 WKB solution↩︎

The amplitude function \(A(\mathcal{Z})\) for each of the nine resistive wave modes is calculated using the ideal amplitude equation (29 ), and the WKB solution is computed by summing the modes as, e.g., \[\label{eqn:uwkb} u_{\text{WKB}} = \Re\left\{C\sum_{n=1}^9 \mathcal{A}_n(\mathcal{Z})u_{0,n}(x, \mathcal{Z})\exp\left(i\theta_n(z) +i k_y y - i t\right)\right\},\tag{33}\] with \[\begin{equation} \mathcal{A}_n(\mathcal{Z}) = c_n A_n(\mathcal{Z}), \end{equation} \begin{equation}\label{eqn:phaseintergral} \theta_n(z) = \int_{z_{b,n}}^{z}k_{z,n}(\mathcal{Z})\mathrm{d}z' + \Delta\theta_n, \end{equation}\tag{34}\] where the lower limit of integration \(z_{b}\) is set to the lower bound of each wave mode (e.g., \(z_{b} = \mathcal{Z}_{t,\text{IGW-0}}/\textit{Fr}\) for IGW-0). Importantly, the amplitude equation (29 ) and eigenvalue problem (32 ) only define \(\mathcal{A}_n\) and \(\theta_n\) up to multiplicative and additive constants, respectively. Thus, relative amplitudes \(c_n\) (which are distinct from the overall complex amplitude \(C\)) and the phase shifts \(\Delta\theta_n\) must be determined independently.

For IGW-0 and IGW-1, \(c_n\) and \(\Delta\theta_n\) are set such that the IGW modes match the forcing in the simulation at \(\mathcal{Z}_0 = 0.225\): \[\begin{equation} \theta_{\text{IGW-0}}(\mathcal{Z}_0/\textit{Fr}) = \theta_{\text{IGW-1}}(\mathcal{Z}_0/\textit{Fr}) -\pi/2 =0, \end{equation} \begin{equation} \mathcal{A}_{\text{IGW-0}}(\mathcal{Z}_0) = \mathcal{A}_{\text{IGW-1}}(\mathcal{Z}_0) = 1. \end{equation}\] We assume total conversion from IGWs to SM waves, and thus set each \(c_n\) for SM-0, SM-1, evan.-0, and evan.-1 such that \[\begin{equation} \mathcal{A}_{\text{SM-0}}(\mathcal{Z}_{t,\text{IGW-0}})= \mathcal{A}_{\text{evan.-0}}(\mathcal{Z}_{t,\text{IGW-0}}) = \mathcal{A}_{\text{IGW-0}}(\mathcal{Z}_{t,\text{IGW-0}}), \end{equation} \begin{equation} \mathcal{A}_{\text{SM-1}}(\mathcal{Z}_{t,\text{IGW-1}})= \mathcal{A}_{\text{evan.-1}}(\mathcal{Z}_{t,\text{IGW-1}}) = \mathcal{A}_{\text{IGW-1}}(\mathcal{Z}_{t,\text{IGW-1}}). \end{equation}\] At the IGW\(\rightarrow\)SM turning points, the phase differences are assumed to match the theoretical values found in the “axisymmetric” case by [33] such that \[\begin{equation} \theta_{\text{SM-0}}(\mathcal{Z}_{t,\text{IGW-0}}/\textit{Fr}) +\pi/2 =\theta_{\text{evan.-0}}(\mathcal{Z}_{t,\text{IGW-0}}/\textit{Fr})= \theta_{\text{IGW-0}}(\mathcal{Z}_{t,\text{IGW-0}}/\textit{Fr}), \end{equation} \begin{equation} \theta_{\text{SM-1}}(\mathcal{Z}_{t,\text{IGW-1}}/\textit{Fr}) +\pi/2 =\theta_{\text{evan.-1}}(\mathcal{Z}_{t,\text{IGW-1}}/\textit{Fr})= \theta_{\text{IGW-1}}(\mathcal{Z}_{t,\text{IGW-1}}/\textit{Fr}) \end{equation}\] Further, we assume total conversion from SM-0 to AW-0, SM-1 to SM-AW-1, and SM-AW-1 to AW-1 with no associated phase shifts: \[\begin{equation} \mathcal{A}_{\text{AW-0}}(\mathcal{Z}_{t,\text{AW-0}})= \mathcal{A}_{\text{SM-0}}(\mathcal{Z}_{t,\text{AW-0}}),\quad\theta_{\text{AW-0}}(\mathcal{Z}_{t,\text{AW-0}}/\textit{Fr})= \theta_{\text{SM-0}}(\mathcal{Z}_{t,\text{AW-0}}/\textit{Fr}), \end{equation} \begin{equation} \mathcal{A}_{\text{SM-AW-1}}(\mathcal{Z}_{AB})= \mathcal{A}_{\text{SM-1}}(\mathcal{Z}_{AB}),\quad \theta_{\text{SM-AW-1}}(\mathcal{Z}_{AB}/\textit{Fr})= \theta_{\text{SM-1}}(\mathcal{Z}_{AB}/\textit{Fr}), \end{equation} \begin{equation} \mathcal{A}_{\text{AW-1}}(\mathcal{Z}_{t,\text{AW-1}})= \mathcal{A}_{\text{SM-AW-1}}(\mathcal{Z}_{t,\text{AW-1}}),\quad \theta_{\text{AW-1}}(\mathcal{Z}_{t,\text{AW-1}}/\textit{Fr})= \theta_{\text{SM-AW-1}}(\mathcal{Z}_{t,\text{AW-1}}/\textit{Fr}). \end{equation}\] Finally, the complex amplitude \(C\) (which encodes the overall amplitude and phase of \(u_{\text{WKB}}\)) is determined by fitting the WKB solution to \(u_{\text{IVP}}\) at \((x,\mathcal{Z}) = (1/3,0.175)\). To construct \(v_{\text{WKB}}\), we simply replace each \(u_{0,n}\) in (33 ) with the corresponding azimuthal velocity eigenfunction \(v_{0,n}\).

References↩︎