Global 3D Gravito-Acoustic Wave Simulations

• Global 3D gravito-acoustic wave simulations are high-fidelity computational models that capture both gravity and acoustic waves in stratified, rotating stellar interiors.
• They utilize the Hybridizable Discontinuous Galerkin method to discretize complex oscillation equations, ensuring high precision and scalability for helioseismic inversion.
• Validated against 2.5D models, these simulations efficiently leverage advanced solvers and compression techniques to model non-axisymmetric heterogeneities in solar structures.

Global 3D gravito-acoustic wave simulations refer to the computational modeling of time-harmonic oscillatory phenomena in strongly stratified, self-gravitating, rotating stellar objects (primarily the Sun), where both gravity waves (g-modes) and acoustic waves (p-modes) are captured within a three-dimensional geometry. These simulations incorporate realistic solar or stellar structure, potentially including photospheric layers and heterogeneous features, and are essential for interpreting seismic data, performing inversions, and advancing the physical understanding of stellar interiors beyond axisymmetric or linearized approximations (Faucher et al., 27 Nov 2025, Pham et al., 26 Aug 2025).

1. Mathematical Formalism and Physical Model

The core of 3D gravito-acoustic wave simulation is the solution of the linearized adiabatic oscillation equations about a hydrostatic, possibly rotating, stellar background. In the Cowling approximation (neglecting perturbation to the gravitational potential), and using a time-harmonic ansatz (typically $e^{-i\omega t}$), the first-order system for displacement $\bxi$ and pressure perturbation $w = \delta p$ is: $\begin{cases} A\,\bxi + \beta_1\,w + \nabla w = \mathbf{g},\ \nabla \cdot \bxi + \beta_2 \cdot \bxi + \varrho w = h, \end{cases}$ where $A$ is a complex matrix differential operator encoding hydrostatic stratification, rotation, gravity, and attenuation; $\beta_1$, $\beta_2$, and $\varrho$ are coefficient fields built from background pressure, density, and sound-speed profiles; and $$Y = \nabla \phi_0 + \pmb{\Omega} \times (\pmb{\Omega} \times x)$$. The vacuum boundary condition on the computational boundary ($\partial B_0$) is enforced as $w + \bxi \cdot \nabla p_0 = 0$.

Buoyancy effects are governed by the profile of the Brunt–Väisälä (buoyancy) frequency,

$$N^2(r) = \phi_0'(r) \left(\alpha_\rho(r) - \frac{\alpha_p(r)}{\Gamma_1}\right),$$

and the operator alternates between elliptic and hyperbolic regions depending on whether $\omega^2 > N^2(r)$ (acoustic-dominated, elliptic) or $\omega^2 < N^2(r)$ (gravity-dominated, hyperbolic).

A first-order "Liouville" reformulation is often used to regularize the surface behavior, mapping

$\bu_L = \sqrt{\rho_0}\, \bxi, \quad w_L = \frac{\delta p}{\sqrt{\rho_0}},$

and leading to a compact mixed system (Faucher et al., 27 Nov 2025).

2. Numerical Methodologies

The prevalent approach for 3D global gravito-acoustic wave simulations employs the Hybridizable Discontinuous Galerkin (HDG) method. HDG discretization partitions the computational domain into tetrahedral elements, with polynomial approximations for both displacement and pressure-like variables local to each element and a set of global hybrid "trace" unknowns on inter-element faces.

Local elemental equations (strong and weak forms) are coupled via interface continuity enforced through numerical fluxes. The stabilization parameter $\tau$, critical for controlling numerical oscillations, is designed based on characteristics of the operator ($\omega^2$ vs. $N^2$) and includes both "convection" and "diffusion" components: $$\tau = |n^\mathrm{T} A^{-1} \beta_1| + \alpha_\mathrm{tune}|n^\mathrm{T} A^{-1} n|, \quad \alpha_\mathrm{tune} = \mathcal{O}(\omega).$$ Static condensation eliminates interior degrees of freedom locally, yielding a much smaller global linear system in the trace unknowns (roughly 20–27% the size of the full DG system) (Faucher et al., 27 Nov 2025).

For the resulting large, sparse, complex-valued linear systems (with up to $\mathcal{O}(10^8)$ unknowns), direct sparse solvers such as MUMPS are used, enhanced by block low-rank (BLR) compression (off-diagonal block SVD truncation), mixed-precision storage (based on singular value spectra), and block-graph symbolic analysis for efficient parallel factorization (Faucher et al., 27 Nov 2025, Pham et al., 26 Aug 2025).

3. Physical Regimes, Boundary Conditions, and Source Modeling

The simulation domain encompasses the full Sun plus its atmosphere, typically modeled up to $r = 1.001 R_\odot$. Boundary conditions are specified for physical fidelity and regularity: the surface is treated as a vacuum, enforcing a Lagrangian pressure-jump condition, and axis or center singularities are excluded through regularity conditions ($\bxi \cdot \mathbf{n} = 0$).

Point sources (modeling localized excitation) are implemented via right-hand-side terms localized to mesh cells, and heterogeneous perturbations (such as active regions or convection snapshots) are encoded as spatial variations in sound speed $c(r, \theta, \phi)$.

In settings with symmetric or axisymmetric background profiles, azimuthal decomposition (Fourier expansion in the $\phi$ direction) enables "2.5D" reductions, where the full 3D system is recast as independent 2D problems for each azimuthal mode $m$—substantially reducing computational cost. For fully asymmetric (heterogeneous) backgrounds, the full 3D problem must be solved (Pham et al., 26 Aug 2025).

4. Validation and Benchmarks

Validation of 3D simulations relies on cross-comparison with 2.5D axisymmetric solutions, particularly for spherically symmetric solar models. When tested at frequencies of 1–2 mHz and for a range of source depths, the relative difference between 3D and 2.5D solutions on meridional planes is better than $10^{-5}$ outside the immediate neighborhood of the Dirac source, confirming the high precision of the 3D HDG approach (Faucher et al., 27 Nov 2025).

Mesh requirements are dictated by the shortest physical wavelength—the highest examined frequencies (e.g., 6 mHz) require $\sim 10$–20 points per wavelength, leading to mesh sizes of up to $1.2 \times 10^6$ tetrahedral cells. BLR and face-block graph analysis streamline direct factorization, reducing memory to about 20% of full-rank requirements and permitting factorizations of hundreds of millions of variables within practical wall times on high-performance computing clusters (Faucher et al., 27 Nov 2025).

5. Applications: Heterogeneities, Nonlinearity, and Helioseismic Inference

Global 3D gravito-acoustic wave solvers permit detailed forward modeling in the presence of physically relevant solar heterogeneities. Experiments imposing wave-speed perturbations from synoptic magnetograms (active regions) show that the Born approximation (linear response) holds up to $\sim$20% amplitude perturbations; for higher amplitudes, nonlinear effects are apparent. Conversely, convection-induced perturbations (amplitudes $\sim 1\%$) create distributed, frequency-dependent modifications to the wavefield but lack isolated signatures (Faucher et al., 27 Nov 2025).

Such simulations are crucial for quantitatively interpreting high-resolution helioseismic data, enabling seismic inversions, imaging of active and far-side regions, and investigation of nonlinear wave-heterogeneity interactions outside the regime of perturbative methods.

6. Related Results in 3D Hydrodynamic Wave Generation Studies

Complementary work in 3D hydrodynamic modeling of massive stars (e.g., Rogers et al. (Edelmann et al., 2019)) focuses on self-consistent convection, IGW excitation, and propagation from convective to radiative zones. In these models, the wavefield is dominated by core-generated plumes, with the kinetic energy spectrum displaying broad-frequency, double power-law behavior and the resulting IGW standing modes matching GYRE eigenvalues to within a few μHz. These simulations also confirm that high-frequency wave breaking seen in 2D is suppressed in 3D due to enhanced damping, and that surface signatures scale as predicted by pseudomomentum conservation, directly linking velocity and temperature fluctuations to observable brightness variations. However, these simulations filter out p-modes via the anelastic approximation (Edelmann et al., 2019).

7. Outlook and Future Directions

Global 3D gravito-acoustic modeling via HDG and advanced direct solvers establishes a scalable, high-fidelity framework for solar and stellar seismology. This approach is extensible to non-axisymmetric and fully heterogeneous backgrounds, and supports the integration of magnetoconvective effects or more advanced non-ideal physical processes at the expense of additional solver and basis complexity (Pham et al., 26 Aug 2025). The methodology directly supports high-resolution seismic inversion, characterization of localized phenomena, and benchmarking of future hybrid or reduced models, positioning it as a crucial tool in the interpretation and quantitative analysis of next-generation helioseismic and asteroseismic observations.