3D Modeling of Solar Oscillations with Hybridizable Discontinuous Galerkin Method

TL;DR

The paper tackles the challenge of forward modeling 3D gravito-acoustic solar oscillations in a non-symmetric, highly stratified background. It advances a time-harmonic, adiabatic formulation discretized with Hybridizable Discontinuous Galerkin (HDG), and introduces Liouville transformations to dramatically improve system conditioning. By combining hp-adaptivity, block low-rank compression, and mixed-precision storage within a MUMPS-based direct solver, the authors achieve scalable, memory-efficient 3D simulations with multiple right-hand sides. Validation against an axisymmetric solver and experiments with 3D wave-speed perturbations (active regions and convection) demonstrate the method’s accuracy and its potential to support non-spherical helioseismic analyses and farside imaging. The work enables high-fidelity forward modeling in 3D with realistic solar backgrounds and paves the way for more accurate interpretation of helioseismic observables in complex solar interiors.

Abstract

With increasing quantity and quality of solar observations, it becomes essential to account for three-dimensional heterogeneities in wave modeling for seismic data interpretation. In this context, we present a 3D solver of the time-harmonic adiabatic stellar oscillation equations without background flows on a domain consisting of the Sun and its photosphere. The background medium consists of 3D heterogeneities on top of a radial strongly-stratified standard solar model. The oscillation equations are solved with the Hybridizable Discontinuous Galerkin (HDG) method, considering a first-order formulation in terms of the vector displacement and the pressure perturbation. This method combines the high-order accuracy and the parallelism of DG methods while yielding smaller linear systems. These are solved with a direct solver, with block low-rank compression and mixed-precision arithmetic to reduce memory footprint. The trade-off between compression and solution accuracy is investigated, and our 3D solver is validated by comparing with resolution under axial symmetry for solar backgrounds. The capacity of the solver is illustrated with wave speed heterogeneities characteristic of two physical phenomena: active regions and convection. We show the importance of global 3D gravito-acoustic wave simulations, in particular when the amplitudes of the perturbations are strong and their effect on the wavefield cannot be estimated by linear approximations.

Figures (17)

• Figure 1: Illustration of the degrees of freedom (dofs) for HDG discretization in two dimensions with polynomial order $\mathfrak{p}=4$. The volume unknowns are represented with the black dofs (unknowns $\boldsymbol{u}_h$ and $w_h$, hence the black dofs have multiplicity $(\mathrm{dimension}+1)$). The green dofs (of multiplicity 1) are for the numerical trace $\lambda$. Only the dofs of the numerical traces are used in the global linear system.
• Figure 2: The standard solar backgrounds from model S are radial christensen1996current, i.e., they only depend on the distance to the origin $r=\Vert\mathbf{x}\Vert$. The models are provided on a scaled interval such that $r=1$ corresponds to the surface. The computations are performed up until the maximal radius ${r_{\max}}=1.001$, Pham2025axisymreport.
• Figure 3: $hp$-adaptivity using $\texttt{mesh06M}$ and frequencies 1.5m (left), 2.0m (middle), and 2.5m (right). The cells are refined near surface, and the polynomial orders varying between 2 and 8 depending on the local-to-the-cell wavelength. For visualization, only the interior edges are shown.
• Figure 4: Size of the HDG global matrix with frequency and mesh (the numbers are given \ref{['table:matrix-analytics']}). The dashed lines correspond to the best quadratic approximation function.
• Figure 5: Simulations using $\texttt{mesh06M}$ at frequency 1m for two different sources. On the left, panels show the solution at a fixed radius $r=1$, and on the right panels are cuts on the $xz$-plane. On the right panels, the top row corresponds to the real part of the solution $w$ and the bottom row to the relative error $\mathfrak{e}_{xz}$ of \ref{['relative-difference:xz']}. The experiments use the radial solar background model \ref{['section:implementation']} and attenuation $\gamma_{\mathrm{att}}/(2\pi)=10.0µ\Hz$.

Theorems & Definitions (2)

• Remark 1
• Remark 2: Model representation on the mesh