Helioseismic Modeling of Background Flows

Since the first measurements of 5 minute solar oscillations nearly 60 yr ago (Leighton et al. 1962; Claverie et al. 1979) helioseismology has evolved into a rich and diverse field using surface observations to probe solar depths, inferring solar parameters and dynamics. Comprehensive overviews of helioseismology can be found in Christensen-Dalsgaard (2002), Di Mauro (2003), and Gizon & Birch (2005). The global nature of these oscillations has allowed for some of the most precise measurements of the solar interior that are currently available, which have become an indispensable tool in measuring internal solar rotation (Duvall et al. 1984; Schou et al. 1998) and meridional circulation (Giles et al. 1997; Zhao & Kosovichev 2004).

The nature of solar rotation is a complex problem and carries with it important information on the dynamic structure of the Sun. Accurate modeling of internal rotation is vital to understanding the effects that rapid rotation rates and angular momentum transport can play in the internal mixing of elements and hence the evolution of solar structure. Observations of the solar surface show a pattern of differential rotation (Snodgrass & Ulrich 1990) which varies widely in its period at the equator (∼24 days) to near the poles (∼30 days). This pattern of differential rotation is mimicked throughout the convective interior (the upper 30% of the solar radius), with a slight maximum in the velocity of subsurface layers (∼0.95R_⊙). These rates begin to converge at the tachocline (Kosovichev 1996) where the differential rotation is coupled to a solid core, rotating at a rate of ∼430 nHz. Even though helioseismology has made unprecedented strides in modeling interior solar dynamics, many unanswered questions remain. Inversion results near the poles show some discrepancies between the Global Oscillation Network Group (GONG) and the Michelson Doppler Imager (MDI; Schou et al. 2002) data pipelines as well as newer HMI (Howe et al. 2011) observations. There are also conflicting measurements of the solar core's rotation rate using low-degree modes (l = 1–4)—ranging from significantly lower rates (BiSON; Chaplin et al. 1996) to much higher ones (IRIS; Lazrek et al. 1996). Recent measurements of g-modes have also implied rotation rates more than twice as fast as previous estimates (Fossat et al. 2017).

Perpendicular to global rotation, we see strong poleward flows (20 m s⁻¹) in each hemisphere (Hathaway 1996). These meridional mass flows operate as mechanisms that distribute angular momentum and magnetic flux throughout the convective interior (Hathaway et al. 1996). Techniques in local helioseismology use perturbations in acoustic travel times to observe these flows (Gizon & Birch 2005), however it becomes more difficult to clearly resolve structures at greater depths. There is still no consensus on the nature and location of the return flow of meridional circulation cells—while commonly thought to sit in the relatively deep region at the base of the tachocline (Giles 1999), new techniques and measurements have begun to question this assumption (Mitra-Kraev & Thompson 2007; Hathaway 2012). The structure of these circulation cells has also been put in doubt, with recent measurements inferring more than just a single cell per hemisphere (Zhao et al. 2013), cf. Gizon et al. (2020).

Through persistent uncertainties on the nature of solar structure, forward modeling offers an important avenue for the validation and systematic analysis of the helioseismology techniques used to infer internal solar parameters. Numerical acoustic solar models have been instrumental in testing these techniques under approximate solar conditions; since the seminal work of Jensen et al. (2003), simulating subsurface sound-speed perturbations (Parchevsky & Kosovichev 2007; Parchevsky et al. 2014), 3D Cartesian models of the solar atmosphere have been employed to validate local time–distance measurements in regimes of convection (Braun et al. 2007) and magnetic fields (Cameron et al. 2008; Khomenko et al. 2009; Parchevsky & Kosovichev 2009a, 2009b; Felipe et al. 2016). The limitations of the computational domains, needed to probe deeper structures, led to the development of global spherical models (Hanasoge et al. 2006, 2007; Hartlep et al. 2007, 2013; Papini et al. 2015; Gizon et al. 2017)—providing key validations of inversion techniques and inferred observations of sound-speed perturbations and background flow structures in the hydrodynamic regime. These models have provided a strong basis for evaluating mean flow structures on the Sun—providing support (Hartlep et al. 2013) for the inference of double-cell (Zhao et al. 2013; Chen 2019) meridional velocity profiles in HMI observations; however, with recent inferences that reaffirm the single-cell structure (Gizon et al. 2020) on MDI and GONG data, it becomes clear that a detailed systematic investigation of these techniques and the limits of resolving internal structures is Required. In order to address some of these issues we demonstrate the efficacy of a new pseudospectral global acoustic algorithm—developed to quickly and flexibly compute stochastically excited oscillations over wide varieties of static or dynamic three-dimensional background velocity fields. Such fully global solar simulations have the potential to offer insights on the global effects of complex hydrodynamic structures in the solar interior—from modeling meridional circulation and center-to-limb effects (Chen 2019) to exploring the information that g-modes carry about the rotating solar core. The computational techniques presented in this paper will also become the basis for adding linear contributions of global magnetic field terms in future work.

This paper is organized as follows. In Section 2 we describe the mathematical background and computational setup of the model. In Section 3 we detail the specific numerical methods we used to compute our governing equations. In Section 4 we reproduce the solar p-mode spectrum and frequency splitting due to differential rotation. We continue our validation by reproducing measurements of a simple single-cell model of meridional circulation using local helioseismology techniques in Section 5. Finally, we discuss future plans for the model and state our concluding remarks in Section 6.

2.1. Mathematical Background

In this section we present the mathematical formulation and computational setup for the Global Acoustic Linearized Euler code—a three-dimensional numerical solver of the linearized compressible Euler equations, used to model acoustic oscillations throughout the solar interior. For the sake of simplicity we assume the adiabatic approximation, where the timescale of heat transfer is much smaller than the period of oscillations and is therefore neglected in the conservation of energy (Christensen-Dalsgaard 2014). The governing equations are solved on a fully global spherical mesh, 0 ≤ ϕ ≤ 2π, 0 ≤ θ ≤ π, 0 ≤ r ≤ R_⊙, in the Cauchy conservation form; they are enumerated as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho ^{\prime} }{\partial t}+{\rm{\Upsilon }}^{\prime} =0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\rm{\Upsilon }}^{\prime} }{\partial t}+{\boldsymbol{\nabla }}:\left({\boldsymbol{m}}^{\prime} \tilde{{\boldsymbol{u}}}+\tilde{\rho }\tilde{{\boldsymbol{u}}}{\boldsymbol{u}}^{\prime} \right)\\ \quad =\,-{{\rm{\nabla }}}^{2}\left(p^{\prime} \right)-{\rm{\nabla }}\cdot \left(\rho ^{\prime} {\tilde{g}}_{r}\hat{{\boldsymbol{r}}}\right)+{\rm{\nabla }}\cdot S\hat{{\boldsymbol{r}}},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial p^{\prime} }{\partial t}=-\displaystyle \frac{{{\rm{\Gamma }}}_{1}\tilde{p}}{\tilde{\rho }}\left(\Space{0ex}{3.25ex}{0ex}{\rm{\nabla }}\cdot \tilde{\rho }{\boldsymbol{u}}^{\prime} +\rho ^{\prime} {\rm{\nabla }}\cdot \tilde{{\boldsymbol{u}}}\right.\\ \qquad -\,\left.\displaystyle \frac{p^{\prime} }{\tilde{p}}\tilde{{\boldsymbol{u}}}\cdot {\rm{\nabla }}\tilde{\rho }+\tilde{\rho }{\boldsymbol{u}}^{\prime} \cdot \displaystyle \frac{{N}^{2}}{g}\hat{{\boldsymbol{r}}}\right).\end{array}\end{eqnarray} \tag{ 3 }$$

In Equations (1)–(3), we solve for oscillations in the potential field (ϕ) independently, and solenoidal contributions are discarded, here ϒ is defined as the divergence of the momentum field m (ϒ = ∇ · m = ∇ · ρ u = ∇² ϕ). ρ, u, p, and g are the density, velocity, pressure, and gravity terms, respectively. Γ₁ is the adiabatic ratio, constant in time in the adiabatic approximation. The full derivation for current form of the energy equation (Equation (3)) can be found in Appendix A. The governing equations are linearized—split into a base flow (tilde) and a perturbation from that base flow (prime), and only the first-order correlation terms are considered; the base values are derived from a theoretical background state, the standard solar model S (Christensen-Dalsgaard et al. 1996). N² is the Brunt–Väisälä frequency where

$$\begin{eqnarray}&&{N}^{2}=g\left(\displaystyle \frac{1}{{{\rm{\Gamma }}}_{1}}\displaystyle \frac{\partial \mathrm{ln}p}{\partial r}-\displaystyle \frac{\partial \mathrm{ln}\rho }{\partial r}\right).\end{eqnarray} \tag{ 4 }$$

Our set of linear equations become unconditionally unstable in regions where this term is negative (N² < 0). These instabilities can be avoided by deriving the conservation of energy (Equation (3)) as a function of the Brunt–Väisälä frequency, allowing us to set slightly negative values in the convection zone to zero directly, without the need to alter background profiles of pressure, density or the adiabatic ratio (Γ₁) as in similar convectively stable models of Hanasoge et al. (2006), Parchevsky & Kosovichev (2007), and Hartlep et al. (2007). Altering the algorithm to maintain stability may introduce small deviations from the original model (Papini et al. 2014), but it is accurate enough to be used in testing helioseismology techniques throughout the convection zone.

The formulation of the source function ($S\hat{{\boldsymbol{r}}}$) is similar to that of Hanasoge et al. (2006), where a thin shallow layer bounded by the model surface (R_⊙) simulates the source of solar oscillations; see Goldreich et al. (1994) and Stein & Nordlund (2001). This function is modeled as a radial Gaussian with a standard deviation of σ = 0.0001R_⊙ ∼ 69.6 km and centered at μ = 0.9995R_⊙—simulating radial force perturbations from the solar surface. The full time-dependent spectrum of the source is generated in frequency space with a Gaussian function centered at μ = 3.2 mHz and with a standard deviation of σ = 1 mHz—simulating the power peak of observed solar oscillations. To mimic a stochastic excitation of our modes (Woodard 1984), we multiply our temporal Gaussian function by a set of random numbers at each frequency interval (f_s = 1/Δt); subsequently applying a Fourier transform to produce a randomized oscillating function for our source, of which a unique one is created for every harmonic degree and azimuthal order (l, m) in our spherical harmonic decomposition.

The material derivative (${\boldsymbol{\nabla }}:\left({\boldsymbol{m}}^{\prime} \tilde{{\boldsymbol{u}}}+\tilde{\rho }\tilde{{\boldsymbol{u}}}{\boldsymbol{u}}^{\prime} \right)$) is solved in its conservation form in order to fully account for effects that background velocities exert on the conservation of momentum (Equation (2)) in the acoustic regime. This is an important requirement for properly simulating the effects of differential rotation and meridional circulation on the frequency of solar oscillations.

2.2. Radial Grid

2.3. Boundary Conditions

The boundary layers are solved as simple reflective walls with a zero velocity perturbation condition (${\boldsymbol{u}}^{\prime} =0$). To avoid nonphysical surface reflections from affecting our solution, we implement a buffer layer by placing a damping factor (σ) into our governing equations:

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}\displaystyle \frac{{\rm{\partial }}{\rho }^{{\rm{{\prime} }}}}{{\rm{\partial }}t} & = & -{{\rm{\Upsilon }}}^{{\rm{{\prime} }}}-\sigma {\rho }^{{\rm{{\prime} }}},\\ \displaystyle \frac{{\rm{\partial }}{{\rm{\Upsilon }}}^{{\rm{{\prime} }}}}{{\rm{\partial }}t} & = & -{{\rm{\nabla }}}^{2}{p}^{{\rm{{\prime} }}}+{ \mathcal Y }-\sigma {{\rm{\Upsilon }}}^{{\rm{{\prime} }}},\\ \displaystyle \frac{{\rm{\partial }}{p}^{{\rm{{\prime} }}}}{{\rm{\partial }}t} & = & { \mathcal P }-\sigma {p}^{{\rm{{\prime} }}}.\end{array}\end{array}\end{eqnarray} \tag{ 5 }$$

To avoid precision errors from our damping term, it is computed implicitly in our time-discretization scheme using the integrating factor method. Damping is initiated at the model surface (R_⊙) and steadily increased into the atmospheric layers, mimicking the escape of acoustic oscillations above the cutoff frequency (>5 mHz).

The time-discretization scheme used in governing Equations (1)–(3) is the second-order accurate (${ \mathcal O }({h}^{2})$) Adams–Bashforth method used to compute the effect of external contributions of background fields to our conservation of momentum and conservation of energy (Equations (2)–(3)) considerations, denoted by ${ \mathcal Y }$ and ${ \mathcal P }$ in Equation (5), respectively. To advance our continuity and conservation of momentum forward in time ($\rho ^{\prime} ({\rm{\Upsilon }}^{\prime} ),{\rm{\Upsilon }}^{\prime} (p^{\prime} )$) (Equations (1)–(2)), we use the implicit first-order accurate (${ \mathcal O }(h)$) backward Euler method to help maintain the stability of our solution.

Spatial differentiation in the radial direction employs the second-order accurate central finite-difference scheme (Equation (6)) when computing the radial component of the Laplacian (∇²) in Equation (2):

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\displaystyle \frac{\partial f}{\partial r}\right)=\displaystyle \frac{2}{{\rm{\Delta }}{r}_{k}+{\rm{\Delta }}{r}_{k+1}}\\ \quad \times \,\left({r}_{k+0.5}^{2}\displaystyle \frac{{f}_{k-1}-{f}_{k}}{{\rm{\Delta }}{r}_{k+1}}-{r}_{k-0.5}^{2}\displaystyle \frac{{f}_{k}-{f}_{k-1}}{{\rm{\Delta }}{r}_{k}}\right),\end{array}\end{eqnarray} \tag{ 6 }$$

where subscript k denotes separate radial mesh points and Δr_k represents the radial distance between points k − 1 and k. The computation of the first radial derivative also employs a second-order accurate polynomial central finite-difference scheme (Figure 7).

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial f}{\partial r} & = & -\displaystyle \frac{{\rm{\Delta }}{r}_{k+1}}{{\rm{\Delta }}{r}_{k}({\rm{\Delta }}{r}_{k+1}+{\rm{\Delta }}{r}_{k})}{f}_{k-1}\\ & & +\,\displaystyle \frac{{\rm{\Delta }}{r}_{k+1}-{\rm{\Delta }}{r}_{k}}{{\rm{\Delta }}{r}_{k+1}{\rm{\Delta }}{r}_{k}}{f}_{k}+\displaystyle \frac{{\rm{\Delta }}{r}_{k}}{{\rm{\Delta }}{r}_{k+1}({\rm{\Delta }}{r}_{k+1}+{\rm{\Delta }}{r}_{k})}{f}_{k+1}.\end{array}\end{eqnarray} \tag{ 7 }$$

The precision of these techniques and their agreement with theoretical calculations of eigenmodes is quite high (Figure 1), serving well to validate the computational techniques used in this algorithm; higher-order accurate schemes may improve our results, but in the large variations in our nonuniform radial regime, many of the improvements that we expect may be greatly diminished.

Zoom In Zoom Out Reset image size

Differentiation tangent to the sphere is done in the pseudospectral regime through spherical harmonic decomposition using the Libsharp spherical harmonic library (Reinecke & Seljebotn 2013). Tangential resolution is defined by the maximum allowable quantum number (l_max(r)) at each radial point. This value controls the mesh size of our Gauss–Legendre grid, containing ${N}_{\phi }=3{l}_{\max }$ azimuthal grid points and ${N}_{\theta }=3{l}_{\max }/2$ latitudinal grid points. These values are chosen to avoid aliasing during spherical harmonic decomposition. The azimuthal mesh points (N_ϕ ) are spaced at even intervals between 0 < ϕ < 2π, while the latitudinal points (N_θ ) are placed at the roots of the corresponding Legendre polynomial between 0 < θ < π. To avoid oversampling at high latitudes the Libsharp library natively implements polar optimization, or the "reduced Gauss–Legendre grid."

To compute the tangential components of the divergence we use vector and tensor spherical harmonic bases (VSH and TSH) to calculate our terms spectrally, as defined in Sections 3.1 and 3.2 respectively.

3.1. Vector Spherical Harmonics

Using the "pure-spin" VSH components defined by Thorne (1980), we can expand an arbitrary vector field (E) into the following linearly independent basis.

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{E}}=\displaystyle \sum _{m=-\mathrm{lmax}}^{\mathrm{lmax}}\displaystyle \sum _{l=| m| }^{\mathrm{lmax}}\left({E}_{{lm}}^{r}(r){{\boldsymbol{Y}}}_{{lm}}\right.\\ \quad \left.+\,{E}_{{lm}}^{(1)}(r){{\boldsymbol{\Psi }}}_{{lm}}+{E}_{{lm}}^{(2)}(r){{\boldsymbol{\Phi }}}_{{lm}}\right),\end{array}\end{eqnarray} \tag{ 8 }$$

where E^r _lm is the radial vector component, and ${E}_{{lm}}^{(1)}$, ${E}_{{lm}}^{(2)}$ are components tangential to the surface of the 2-sphere. Transformations between the spherical coordinate basis and our VSH basis is achieved by a set of functions defined in Novak et al. (2010), and computed using recurrence relations of the spherical harmonic (Y_lm ). We define the divergence of our vector field ( E in Equation (8)) in our new basis with Equation (9):

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{E}}=\displaystyle \sum _{m=-\mathrm{lmax}}^{\mathrm{lmax}}\displaystyle \sum _{l=| m| }^{\mathrm{lmax}}\left(\displaystyle \frac{{{dE}}_{{lm}}^{r}}{{dr}}+\displaystyle \frac{2}{r}{E}_{{lm}}^{r}-\displaystyle \frac{l(l+1)}{r}{E}_{{lm}}^{(1)}\right){Y}_{{lm}}.\end{eqnarray} \tag{ 9 }$$

The radial derivative is computed using the finite-difference method described in Section 3.

3.2. Tensor Spherical Harmonics

The TSH basis is built on groups 0 and 2 of the irreducible representation of rotation in SO(3) (${{ \mathcal D }}^{0}$, ${{ \mathcal D }}^{2}$ where $t={{ \mathcal D }}^{0}\,+{{ \mathcal D }}^{1}+{{ \mathcal D }}^{2}$), which correspond to the trace and the symmetric traceless tensor respectively (Mathews 1962). This basis is coupled with spin-0 and spin-2 spherical harmonics to form six irreducible representations of a symmetric tensor, which can be rearranged into the "pure-spin tensor harmonics" defined by Zerilli (1970) and Thorne (1980). We can use these components to expand an arbitrary symmetric tensor ( t ) into the following orthogonal basis (Novak et al. 2010):

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{t}}=\displaystyle \sum _{m=-\mathrm{lmax}}^{\mathrm{lmax}}\displaystyle \sum _{l=| m| }^{\mathrm{lmax}}\left({L}_{{lm}}^{0}{{\boldsymbol{T}}}_{{lm}}^{{{\boldsymbol{L}}}_{0}}+{E}_{{lm}}^{1}{{\boldsymbol{T}}}_{{lm}}^{{{\boldsymbol{E}}}_{1}}\right.\\ \ \ \left.+\,{B}_{{lm}}^{1}{{\boldsymbol{T}}}_{{lm}}^{{{\boldsymbol{B}}}_{1}}+{T}_{{lm}}^{0}{{\boldsymbol{T}}}_{{lm}}^{{{\boldsymbol{T}}}_{0}}+{E}_{{lm}}^{2}{{\boldsymbol{T}}}_{{lm}}^{{{\boldsymbol{E}}}_{2}}+{B}_{{lm}}^{2}{{\boldsymbol{T}}}_{{lm}}^{{{\boldsymbol{B}}}_{2}}\right),\end{array}\end{eqnarray} \tag{ 10 }$$

where L⁰ _lm is the fully radial component (t^rr ) and T⁰ _lm is the transverse (t^{θ θ} , t^{ϕ ϕ} ) portion of the trace. E¹ _lm and B¹ _lm represent the mixed radial/transverse components (t^{r θ} , t^{r θ} ), and E² _lm and B² _lm are symmetric transverse traceless components. We can use this basis to solve for the divergence of a tensor directly into VSH coordinates, as shown in Equation (11).

$$\begin{eqnarray}{\rm{\nabla }}\cdot {\boldsymbol{t}}=\left\{\begin{array}{ll}{\left[{\rm{\nabla }}\cdot {\boldsymbol{t}}\right]}_{{lm}}^{r}{Y}_{{lm}}= & \displaystyle \sum _{m=-\mathrm{lmax}}^{\mathrm{lmax}}\displaystyle \sum _{l=| m| }^{\mathrm{lmax}}\left[\tfrac{1}{{r}^{2}}\tfrac{\partial }{\partial r}\left({r}^{2}{L}_{{lm}}^{0}\right)-\tfrac{1}{r}\left(l(l+1){E}_{{lm}}^{1}+{T}_{{lm}}^{0}\right)\right]{Y}_{{lm}},\\ {\left[{\rm{\nabla }}\cdot {\boldsymbol{t}}\right]}_{{lm}}^{(1)}{Y}_{{lm}}= & \displaystyle \sum _{m=-\mathrm{lmax}}^{\mathrm{lmax}}\displaystyle \sum _{l=| m| }^{\mathrm{lmax}}-l(l+1)\left[\tfrac{1}{{r}^{3}}\tfrac{\partial }{\partial r}\left({r}^{3}{E}_{{lm}}^{1}\right)+\tfrac{1}{r}\left(\tfrac{{T}_{{lm}}^{0}}{2}-(l-1)(l+2){E}_{{lm}}^{2}\right)\right]{Y}_{{lm}},\\ {\left[{\rm{\nabla }}\cdot {\boldsymbol{t}}\right]}_{{lm}}^{(2)}{Y}_{{lm}}= & \displaystyle \sum _{m=-\mathrm{lmax}}^{\mathrm{lmax}}\displaystyle \sum _{l=| m| }^{\mathrm{lmax}}-l(l+1)\left[\tfrac{1}{{r}^{3}}\tfrac{\partial }{\partial r}\left({r}^{3}{B}_{{lm}}^{1}\right)-\displaystyle \frac{(l-1)(l+2)}{r}{B}_{{lm}}^{2}\right]{Y}_{{lm}}.\end{array}\right.\end{eqnarray} \tag{ 11 }$$

Coordinate transformations between the spherical, VSH, and TSH bases is done using relations defined in Novak et al. (2010). We can plug the results of this computation into Equation (9) to compute the divergence dyad (∇: ) in our material derivative (Equation (2)).

In order to validate the mathematical setup (Section 2) and computational techniques (Section 3) described in this paper, we start by matching the generated power spectrum from our model with a theoretical prediction computed using the standard solar model S (Section 4.1). We test the computation of our material derivative (Equation (2)) in two simple regimes of background velocity flows: differential rotation (Section 4.2) and a single-cell model of meridional circulation (Section 5).

The parameters used in this validation are as follows: the surface resolution of the model is set by a spherical harmonic degree of ${l}_{\max }=200$, corresponding to N_ϕ = 600 longitudinal and N_θ = 450 latitudinal mesh points. These values are chosen in order to fully resolve the acoustic modes that intersect the base of the tachocline (a depth of ∼0.70R_⊙), letting us infer flows from any point in the convective interior using local helioseismology techniques (Gizon & Birch 2005). Our simulation is run for a period of 65 hr model time. This is too short to properly resolve flow velocities on the order of those seen in meridional circulation (∼14 m s⁻¹). In order to achieve a desirable signal-to-noise ratio (S/N) in our validation of meridional circulation measurements, we increase the magnitude of our meridional velocity by a factor of 36, to a maximum of 500 m s⁻¹—simulating a S/N that would be observed over ∼9.5 yr, approaching the decade minimum that is estimated to be needed in order to resolve deep meridional flows (Braun & Birch 2008).

4.1. Power Spectrum

Acoustic oscillations on the solar surface can be decomposed into eigenmodes of radial velocity (${u}_{r}^{{\prime} }\hat{{\boldsymbol{r}}}$), representing standing waves throughout the solar interior. These modes can be conceptualized as a combination of the scalar spherical harmonic (Y_lm ) with a frequency dependent radial order (ξ_n ),

$$\begin{eqnarray}&&{u}_{r}^{{\prime} }(r,\theta ,\phi ,t)=\displaystyle \sum _{n,l}\displaystyle \sum _{m=-l}^{l}{\xi }_{n}(r){Y}_{{lm}}(\theta ,\phi ){e}^{i\omega t}.\end{eqnarray} \tag{ 12 }$$

The power spectrum of these modes can be visualized using an l−ν diagram, showing continuous radial modes throughout the solar interior as a function of their frequency (ω = 2πν) and spherical harmonic degree (l) (see Figure 1). The eigenmodes are excited in the frequency range determined by our source function (S, see Section 2.1). By setting negative values of the Brunt–Väisälä frequency to zero we remove the convective instabilities that act as sources of acoustic perturbations normally seen below 2 mHz. We see an agreement in the structure of the eigenmodes generated by our model to theoretical calculations from the model S (Christensen-Dalsgaard et al. 1996), denoted by by the dashed blue lines in Figure 1.

4.2. Rotation

In order to test the computation of our material derivative (Equation (2)), we implement a simple model of differential rotation in a nonrotating reference frame by defining our background velocity term as the angular frequency (Ω) derived from the mean-field model M1 described in Pipin & Kosovichev (2019) and shown in Figure 2.

Zoom In Zoom Out Reset image size

This simple azimuthal velocity flow field model creates easily detectable rotational splitting in the structure of our eigenmodes. The rotational profile will shift up the frequency of prograde modes and shift down the frequency of retrograde modes as a function of their azimuthal order (m). In order to visualize this shift we can use an m−ν diagram of the power spectrum for spherical harmonic degree l = 180 (Figure 3), where our simulated modes reproduce the characteristic tilt due to the average rotation along with the curvature created by the differential rotation in the convection zone. The blue dashed lines show frequency splittings calculated using heliosesimic sensitivity kernels (Schou et al. 1998) for the M1 model (Pipin & Kosovichev 2019).

Zoom In Zoom Out Reset image size

The computation of our material derivative (Equation (2)) can be tested in a regime of global meridional flows by using a simple single-cell model of meridional circulation as our background velocity term ($\tilde{{\boldsymbol{u}}}$, Figure 5). We chose the model described and tested by Hartlep et al. (2013), recreating their measurements as a simple validation of our computational techniques.

In order to infer flows in the model interior using surface measurements, we apply a method of deep focusing (Zhao et al. 2009)—a local helioseismology technique (for an in-depth review, see Gizon & Birch 2005) which uses travel times of acoustic waves to probe structures in solar and stellar interiors. This method consists of choosing two points on the model surface, separated by some angular distance (Δ); the signals at these points are cross-correlated, measuring acoustic travel times of internal oscillations (p-modes). As these waves travel through the resonant cavity of the convective interior they are advected by internal mass flows. Sampling waves traveling in opposite directions results in travel-time differences (δτ) which provide a basis for inferring background velocities. The ray-path approximation under the assumption of Fermat's principle offers a simple relation between these values (Giles 1999).

$$\begin{eqnarray}&&\delta \tau =-2{\int }_{{{\rm{\Gamma }}}_{0}}\displaystyle \frac{\tilde{{\boldsymbol{u}}}\cdot {\boldsymbol{n}}}{{c}^{2}}{ds}.\end{eqnarray} \tag{ 13 }$$

Travel-time differences can be estimated along the unperturbed ray path (Γ₀).

We employ this technique to infer meridional velocities by measuring travel-time differences between southward and northward traveling waves (δ τ_NS). Using each pixel in our data set as a center point, we remap the surrounding 60° × 60° patch into azimuthal equidistant coordinates (Postel's projection) using cubic Hermite splines. The remapped resolution is approximately 0.6° per pixel, the same spatial resolution as our Gauss–Legendre grid (Section 3). An illustration of our method can be seen in Figure 4.

Zoom In Zoom Out Reset image size

A series of six concentric great circles are drawn at diameters of Δ = 12°, 18°, 24°, 30°, 36° and 42°. Sets of pixels (two pixels in width) are selected along the rings in 30° wide northern and southern sectors. The pixels in each sector are averaged together and the signals in opposing sectors are cross-correlated. This process is repeated for every grid point in our model (N_ϕ , N_θ ) and the cross-correlated signal is averaged over every point in the longitude (N_ϕ ) and ±3° (±5 pixels in N_θ ) in the latitude. To further smooth our data, we average the diameter of each great circle over ±2.4°—travel times of the varying diameter distances are interpolated to an estimated time offset based on ray-path theory (Giles 1999). The radial turning points of acoustic waves corresponding to each angular distance are ∼0.93, 0.89, 0.85, 0.81, 0.77, and 0.72 R_⊙, respectively, allowing us to probe the entirety of the convective interior. In Figure 5 we show the meridional flow profile used in the model with dashed lines in the upper hemisphere corresponding to the ray paths for each great circle.

Zoom In Zoom Out Reset image size

5.1. Results

The travel-time differences (δ τ_NS) for each ring diameter are plotted as a function of latitude in Figure 6. These values are compared to theoretical travel-time differences (dashed lines in Figure 6) computed using the ray-path approximation (Equation (13)) employing the standard solar model S (Christensen-Dalsgaard et al. 1996).

Zoom In Zoom Out Reset image size

The travel-time differences for all ring diameters show solid agreement with theoretical predictions as well as the analysis of Hartlep et al. (2013). These results show a key validation of the numerical procedure used to compute the model as well as the deep-focusing techniques used to analyze the data. The error (σ_NS) is calculated using a separate model with no background flows; this reference model uses an identical source function (S, Equation (2)) and analysis sequence, producing the same error profile we see for each ring diameter (Figure 6). We characterize this error as the standard deviation of travel-time differences (δτ) from zero in our reference model, taking the rms over latitudinal grid points.

$$\begin{eqnarray}&&{\sigma }_{\mathrm{NS}}=\sqrt{\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}\delta {\tau }_{i}^{2}}.\end{eqnarray} \tag{ 14 }$$

The characteristic noise profiles seen in models without flows can be subtracted from measured travel-time differences (Figure 6) in order to remove some of the most significant impacts of realization noise on our measurements (Hanasoge et al. 2007). This method of attenuating noise can provide us an opportunity to compare measurements made by computational helioseismology to estimates of the ray-path approximation. The resulting N–S travel-time profile can be found in Appendix B, Figure 9.

In order to further increase the S/N in our measurements we apply a phase-speed filter, defined by a Gaussian function with a width of σ = 0.05v_p , where v_p = ω/L is the phase speed (Nigam et al. 2007). After the application of this filter, the travel-time differences display high precision but do show a latitude-independent systematic offset for different ring diameters. This offset, which seems to have a maximum of approximately ±1 s (Figure 7), is a characteristic of the realization noise in our model and can be removed by subtracting the phase-filtered noise profile for the same model with no background flows. The resulting profiles, which can be found in Appendix B, Figure 10, show remarkable agreement with predictions made using the ray-path approximation (Equation (13)), implying that nonlinear effects do not seem to have dramatic impacts on our measurements of travel times.

Zoom In Zoom Out Reset image size

The error function (σ_NS, Equation (14)) provides a solid foundation for the characterization of realization noise at various acoustic travel depths. We show the error as a function of travel distance (Δ) in Figure 8, along with travel-time differences (δτ) from our reference model without background flows for five separate latitudinal averages (30°N–50°N, 10°N–30°N, 10°S–10°N, 10°S–30°S, and 30°S–50°S) and compare them with the error in the unfiltered signal (Figure 6). The cause of the apparent systematic error is unclear, however, the latitude-dependent offset profiles are strongly linked to the structure of our source (S, Equation (2)), with different random number seeds generating different error profiles. This effect deserves its own systematic investigation for varying parameters of source locations and structures.

Zoom In Zoom Out Reset image size

In Figure 8, we see a similar systematic error structure in both the filtered and unfiltered signals. As we move toward greater depths, however, noise in the unfiltered signal grows significantly, concealing any potential offset. These results may have interesting implications for measuring meridional flow structures at the base of the tachocline. The application of our phase-speed filter seems to have preserved signal quality relatively evenly throughout the convection zone, offering encouraging results for probing flows deep in the solar interior.

We present a global linearized acoustic algorithm with new computational methods that will facilitate the testing and validation of local and global helioseismology techniques in diverse regimes of three-dimensional flows. While helioseismology has been an indispensable tool in exploring the interior dynamics of the Sun, it can have trouble resolving exact profiles of flow, especially at greater depths. Forward modeling offers the opportunity to test the impact of subtle differences generated by a variety of theoretical models of mean mass flows, forming a basis to better interpret observational oscillation data.

The model is validated for two distinct profiles: differential rotation and meridional circulation. These two regimes are critical for understanding and simulating the distribution of angular momentum and magnetic flux that governs the solar magnetic cycle. To simulate these structures, we present an application of pseudospectral techniques which will be necessary components for future models to efficiently compute spherical harmonic resolutions of ${l}_{\max }\gt 300$, previously considered to be too computationally expensive. These resolutions will be necessary to model local helioseismology techniques on sound-speed perturbations due to small-scale structures on the solar surface, creating a link to linear effects of global perturbations.

In future work, we plan to use this model to test differences in acoustic travel-time signatures between models of one and two-cell meridional circulation. The current state surrounding the nature of meridional velocity profiles remains uncertain, with similar techniques producing widely varying results between observational data from HMI (Zhao et al. 2013) and MDI/GONG (Gizon et al. 2020). While both studies employed validation through the use of forward modeling, the theoretical profiles of meridional circulation structures for one and two-cell models varied extremely. A new efficient and flexible code presents a timely opportunity to run a differential study on much more subtle differences between one and two-cell circulation models for both theoretical profiles as well as ones derived from solar simulations. Our investigation (S. Kosovichev 2021, in preparation) will explore subtle differences in flow velocities near the base of the tachocline, focusing on the limits what can be resolved within the two-decade-long window of current solar observations, with and without preparing data using phase-speed filtering. These models will be used for an analysis of realization noise in simulations of stochastic excitations of surface oscillations for large numbers of source profiles and measurement times. This level of systematic investigation will provide a baseline for what can be interpreted from observational data. In further work, we also plan to create a new generation of models that include the linear effects of magnetic field structures on acoustic oscillations. The numerical method demonstrated in this paper will be the basis for computing this effect.

A.M.S. would like to thank the heliophysics modeling and simulation team at NASA Ames Research Center for very productive meetings. This work is supported by the NASA grants 80NSSC19K0630, 80NSSC19K1436, NNX14AB70G, 80NSSC20K0602, and NNX17AE76A.

In our formulation, we employ the adiabatic approximation, where we assume that heat transfer is negligible in the timescale of acoustic oscillations. The relationship between pressure and density in this regime (see Christensen-Dalsgaard 2014) can be expressed in Eulerian form as

$$\begin{eqnarray}&&\displaystyle \frac{\partial p}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}p=\displaystyle \frac{{{\rm{\Gamma }}}_{1}p}{\rho }\left(\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}\rho \right).\end{eqnarray} \tag{ A1 }$$

Assuming that the adiabatic ratio (Γ₁) remains constant, we can linearize Equation (A1) as follows, where primes denote a perturbation from a base (tilde) value.

$$\begin{eqnarray}&&\begin{array}{l}\tilde{\rho }\displaystyle \frac{\partial p^{\prime} }{\partial t}+\rho ^{\prime} \displaystyle \frac{\partial \tilde{p}}{\partial t}+\left(\rho ^{\prime} \tilde{{\boldsymbol{u}}}\cdot {\rm{\nabla }}\tilde{p}+\tilde{\rho }\left(\tilde{{\boldsymbol{u}}}\cdot {\rm{\nabla }}p^{\prime} +{\boldsymbol{u}}^{\prime} \cdot {\rm{\nabla }}\tilde{p}\right)\right)\\ \ \ =\,{{\rm{\Gamma }}}_{1}\left(\tilde{p}\displaystyle \frac{\partial \rho ^{\prime} }{\partial t}+p^{\prime} \displaystyle \frac{\partial \tilde{\rho }}{\partial t}+\left(p^{\prime} \tilde{{\boldsymbol{u}}}\cdot {\rm{\nabla }}\tilde{\rho }+\tilde{p}\left(\tilde{{\boldsymbol{u}}}\cdot {\rm{\nabla }}\rho ^{\prime} +{\boldsymbol{u}}^{\prime} \cdot {\rm{\nabla }}\tilde{\rho }\right)\right)\right).\end{array}\end{eqnarray} \tag{ A2 }$$

In order to avoid convective instabilities in our model we can rewrite Equation (A2) in terms of the Brunt–Väisälä frequency (N², Equation (4)). By plugging in continuity (Equation (1)) and rearranging terms, we are left with:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{\tilde{p}{{\rm{\Gamma }}}_{1}}\displaystyle \frac{\partial p^{\prime} }{\partial t}+\displaystyle \frac{1}{\tilde{\rho }}{\rm{\nabla }}\cdot \left(\tilde{\rho }{\boldsymbol{u}}^{\prime} +\rho ^{\prime} \tilde{{\boldsymbol{u}}}\right)=-{\boldsymbol{u}}^{\prime} \\ \cdot \left(\displaystyle \frac{1}{{{\rm{\Gamma }}}_{1}\tilde{p}}{\rm{\nabla }}\tilde{p}-\displaystyle \frac{1}{\tilde{\rho }}{\rm{\nabla }}\tilde{\rho }\right)-\tilde{{\boldsymbol{u}}}\cdot \left(-\displaystyle \frac{1}{\tilde{\rho }}\left({\rm{\nabla }}\rho ^{\prime} +\displaystyle \frac{p^{\prime} }{\tilde{p}}{\rm{\nabla }}\tilde{\rho }\right)\right).\end{array}\end{eqnarray} \tag{ A3 }$$

Substituting (N², Equation (4)) and rearranging terms further leaves us with the final form of our governing relation for pressure (Equation (3)).

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial p^{\prime} }{\partial t}=-\displaystyle \frac{{{\rm{\Gamma }}}_{1}\tilde{p}}{\tilde{\rho }}\left(\Space{0ex}{3.25ex}{0ex}{\rm{\nabla }}\cdot \tilde{\rho }{\boldsymbol{u}}^{\prime} \right.\\ \quad \left.+\,\rho ^{\prime} {\rm{\nabla }}\cdot \tilde{{\boldsymbol{u}}}-\displaystyle \frac{p^{\prime} }{\tilde{p}}\tilde{{\boldsymbol{u}}}\cdot {\rm{\nabla }}\tilde{\rho }+\tilde{\rho }{\boldsymbol{u}}^{\prime} \cdot \displaystyle \frac{{N}^{2}}{g}\hat{{\boldsymbol{r}}}\right).\end{array}\end{eqnarray} \tag{ A4 }$$

In this section, we show travel-time difference measurements with their noise profiles subtracted (see Section 5.1). Figure 9 shows denoising performed on results presented in Figure 6, while Figure 10 shows the denoising for results presented in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size