Core convection and internal gravity waves in massive stars

Philipp Edelmann

collaborators:
Tami Rogers, Rathish Ratnasingam (Newcastle)
Conny Aerts, Dominic Bowman, May Gade Pedersen (KU Leuven)
Fritz Röpke, Leo Horst (HITS Heidelberg)

Internal gravity waves (IGW) in intermediate-mass/massive stars

  • main sequence: convective core, radiative envelope, (subsurface convection)
  • IGWs excited at convective–radiative interface
  • transport energy, angular momentum, chemical composition
  • coherent modes in cavity
  • waves are amplified towards surface
  • wave breaking at sufficient amplitudes

Rogers+ (2013)

Asteroseismology

  • … probes interior of stars through the patterns that waves leave on the surface
  • … analyzes time-series data of a point-like object
  • modeling by applying oscillation codes to stellar evolution models
  • plenty of observational data from CoRoT, Kepler, TESS
  • … mostly relies on 1D stellar models
  • … relies on assumptions about excitation of waves at the core boundary
  • … nonlinear wave interactions
asteroseismology

source: NASA

The Low-Mach Problem

  • Mach number: $Ma = \frac{u}{c}$
  • stellar interiors are usually at low Mach numbers
  • speed of sound $c = \sqrt{\gamma \frac{p}{\rho}} \propto \sqrt{\frac{T}{\mu}}$ (ideal gas)
  • many Godunov-type schemes show problems for $Ma \rightarrow 0$
  • solution dominated by numerical dissipation (e.g. Guillard & Viozat, 1999)
Solutions
  • control numerical viscosity for $Ma \rightarrow 0$
  • modify the PDEs to remove problematic terms
    (various sound-proof approaches)
  • use a different type of discretization (possibly more viscous)

Gresho Vortex

Gresho & Chan (1990); Liska & Wendroff (2003)

  • rotating 2D vortex, stabilized by pressure gradient
  • stationary solution to the Euler equations
  • Mach number adjustable by $p_0$
$$ u_{\phi} (r) = \begin{cases} 5 r & 0 \le r < 0.2 \\ 2-5 r & 0.2 \le r < 0.4 \\ 0 & 0.4 \le r \end{cases} $$
$$ p(r) = \begin{cases} p_0 + \frac{25}{2} r^2 & 0 \le r < 0.2 \\ p_0 + \frac{25}{2} r^2 + 4 (1-5 r-\ln 0.2 + \ln r) & 0.2 \le r < 0.4 \\ p_0 - 2 + 4 \ln 2 & 0.4 \le r \end{cases} $$

Gresho Vortex

Standard Roe Scheme

Preconditioned Roe Scheme

(Miczek+, 2015)

$$ \vec{F}_{i+1/2} = \frac{1}{2} \left( \vec{F}(\vec{U}^L_{i+1/2}) + \vec{F}(\vec{U}^R_{i+1/2}) - (\color{red}{P^{-1}}|\color{red}{P}A|)_\text{roe}(\vec{U}^R_{i+1/2} - \vec{U}^L_{i+1/2}) \right) $$ $$ P_\vec{V} = \begin{pmatrix} 1 & n_x \frac{\rho\delta\Mr}c{} & n_y \frac{\rho\delta\Mr}c{} & n_z \frac{\rho\delta\Mr}c{} & 0 & 0\\ 0 & 1 & 0 & 0 & -n_x\frac{\delta}{\rho c \Mr} & 0\\ 0 & 0 & 1 & 0 & -n_y\frac{\delta}{\rho c \Mr} & 0\\ 0 & 0 & 0 & 1 & -n_z\frac{\delta}{\rho c \Mr} & 0\\ 0 & n_x \rho c \delta \Mr & n_y \rho c \delta \Mr & n_z \rho c \delta \Mr & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} $$ $$\delta = \frac{1}{\min(1, \max(M,M_\text{cut}))} - 1$$

Gresho vortex

Barsukow+ (2017)

Gresho vortex

kinetic energy

Seven-League Hydro code

  • solves the compressible Euler equations in 1-, 2-, 3-D
  • explicit and implicit time integration
  • flux preconditioning to ensure correct behavior at low Mach numbers
  • other low Mach number schemes (e.g. AUSM$^+$-up)
  • works for low and high Mach numbers on the same grid
  • hybrid (MPI, OpenMP) parallelization (works up to 458 752 cores)
  • several solvers for the linear system:
    BiCGSTAB, GMRES, Multigrid, (direct)
  • arbitrary curvilinear meshes
    using a rectangular computational mesh
  • gravity solver (monopole, Multigrid)
  • radiation in the diffusion limit
  • general equation of state
  • general nuclear reaction network

Planar gravity wave test

Horst+ (in prep.)

The anelastic approach

Our 3D code @ Newcastle
  • working title: SPIN (Stellar and Planetary INteriors)
  • anelastic: small deviations from hydrostatic background, no sound waves
  • optional MHD
  • explicit viscosity, thermal diffusivity, resistivity
  • spherical harmonics discretisation in angular direction
  • pseudo-spectral (nonlinear terms computed in real space)
  • finite differences in radial direction
  • part implicit (linear terms), part explicit (nonlinear terms)
  • ideal gas, not tracing composition (yet)
  • similar to ASH, Rayleigh
  • Glatzmaier (2013) is an excellent guide to these methods

Some properties of the anelastic simulations

  • high diffusivity
  • higher luminosity
  • anelastic approximation
  • decomposition into spherical harmonics
  • hard to cover many $\rho$ scale heights
  • scaling to many cores is hard
available resolution forces us to use $\kappa$ and $\nu$ much higher than the stellar values
compensate high diffusivity with higher luminosity
$$v_\text{conv} \propto L^{1/3}$$
\begin{align} \renewcommand{\vec}{\mathbf} \newcommand{\PD}[2]{\frac{\partial #1}{\partial #2}} \label{eq:anelastic-rho} \nabla \cdot \overline{\rho} \vec{v} = 0, \end{align} \begin{align} \label{eq:anelastic-v} \PD{\vec{v}}{t}&= - (\vec{v} \cdot \nabla) \vec{v} %-\nabla P - \left(T - \overline{\left(\PD{T}{p}\right)}p\right) \frac{\overline{g}}{\overline{T}}\vec{\hat{r}}\\ - \nabla P - C \overline{g} \vec{\hat{r}} + 2(\vec{v} \times \vec{\hat{z}} \Omega) \\ %- v_r \PD{\Omega(r)}{r} \hat{\phi} \nonumber & + \overline{\nu} \left( \nabla^2 \vec{v} + \frac{1}{3} \nabla (\nabla \cdot \vec{v}) \right),\\ \label{eq:anelastic-T} \PD{T}{t}&= - (\vec{v} \cdot \nabla) T + (\gamma - 1) T h_\rho v_r\\ \nonumber %& - \gamma \left[ \frac{d\overline{T}}{dr} - \left(\frac{d\overline{T}}{dr}\right)_\text{ad}\right] v_r\\ &- v_r \left( \PD{\overline{T}}{r} - (\gamma - 1) \overline{T} h_\rho \right) + \frac{\overline{Q}}{c_v \overline{\rho}}\\ \nonumber & + \frac{1}{c_v\overline{\rho}} \nabla \cdot (c_p \overline{\kappa}\overline{\rho}\nabla T) + \frac{1}{c_v\overline{\rho}} \nabla \cdot (c_p \kappa_r\overline{\rho}\nabla \overline{T}).%\\ %+ (\gamma - 1) T h_\rho v_r + \gamma \overline{\kappa} \left(\nabla^2 T + h_\rho \PD{T}{r}\right)\\ %\nonumber %& %+ \gamma \overline{\kappa} \left(\nabla^2 T + h_\rho \PD{\overline{T}}{r}\right) + \frac{\overline{Q}}{c_v}. \end{align}

no sound waves (or p modes) possible in simulation
can easily extract different $(l,m)$ components of the variables
simulation domain:
$0.01\,R_\star$ to $0.9\,R_\star$

strong scaling on Pleiades at NASA Ames

Our star

Edelmann+ (2019)

  • MESA model
  • $3\,\mathrm{M}_\odot$ ZAMS (B-type)
  • $Z=10^{-2}$
  • 1500 radial zones (400 in CZ)
  • angular resolution:
    3741 modes ($128 (\theta) \times 256 (\phi)$) or
    14706 modes ($256 (\theta) \times 512 (\phi)$)
Wave signatures in envelope
What do the observations look like

Bowman et al. (2019)

SLH simulations (Leo Horst, HITS Heidelberg)
  • 2D so far (equatorial annulus)
  • no explicit viscosity, stellar thermal diffusivity
  • fully compressible
  • $L=10^3\,L_\star$

Horst+ (in prep.)

Code-comparison project

  • started at Stellar Hydro Days V (2019)
  • simple convection problem in a Cartesian box, ideal gas
  • well-defined outputs: power spectrum, concentration of mixed fluids, …
  • two setups: just convection zone, convective–radiative interface
  • participating codes so far: SLH, PPMstar, MAESTROeX, Castro, Dedalus, Pencil, …
  • more information: https://sites.google.com/view/stellarhydrodays
    or: philipp@slh-code.org

credit: Robert Andrassy (HITS)

Some advertising

Conclusions
  • low Mach numbers pose (solvable) numerical challenges in accuracy and efficiency
  • IGW propagation can be subtly changed by the numerical scheme
  • convection excites modes and stochastic IGWs that travel to the surface of massive stars
  • observations show low-frequency signal with characteristic slope that could be of IGW origin
Outlook
  • lower viscosity/diffusivity $\rightarrow$ closer to stellar $L$, lower $f$
  • simulation of radiation zone only (impose waves through inner boundary)
  • effect of $B$ field on core convection,wave excitation, propagation
Rayleigh code
  • 3D pseudo-spectral code
  • 2D domain decomposition
  • original developer: Nick Featherstone (CU Boulder)
  • efficient scaling up to 10⁴ cores (Matsui+, 2016)
  • GPLv3 licenced
  • github.com/geodynamics/Rayleigh
  • custom reference states (e.g. from MESA)
  • Cartesian and fully compressible versions planned