Fully compressible simulations of waves and core convection in main-sequence stars
L. Horst, P. V. F. Edelmann, R. Andrassy, F. K. Roepke, D. M. Bowman, C. Aerts, R. P. Ratnasingam
AAstronomy & Astrophysics manuscript no. oscillations c (cid:13)
ESO 2020June 5, 2020
Fully compressible simulations of waves and core convectionin main-sequence stars
L. Horst , P. V. F. Edelmann , , R. Andrássy , F. K. Röpke , , D. M. Bowman , C. Aerts , , , and R. P. Ratnasingam Heidelberger Institut für Theoretische Studien, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germanye-mail: [email protected] School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne, NE1 7RU, UK X Computational Physics (XCP) Division and Center for Theoretical Astrophysics (CTA), Los Alamos National Laboratory, LosAlamos, NM 87545, USA Zentrum für Astronomie der Universität Heidelberg, Institut für Theoretische Astrophysik, Philosophenweg 12, 69120 Heidelberg,Germany Instituut voor Sterrenkunde, KU Leuven, Celestijnenlaan 200D, 3001, Leuven, Belgium Department of Astrophysics, IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL, Nijmegen, The Netherlands Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, GermanyReceived ; accepted
ABSTRACT
Context.
Recent, nonlinear simulations of wave generation and propagation in full-star models have been carried out in the anelasticapproximation using spectral methods. Although it makes long time steps possible, this approach excludes the physics of sound wavescompletely and rather high artificial viscosity and thermal di ff usivity are needed for numerical stability. Direct comparison with ob-servations is thus limited. Aims.
We explore the capabilities of our compressible multidimensional hydrodynamics code SLH to simulate stellar oscillations.
Methods.
We compare some fundamental properties of internal gravity and pressure waves in 2D SLH simulations to linear wavetheory using two test cases: (1) an interval gravity wave packet in the Boussinesq limit and (2) a realistic 3 M (cid:12) stellar model with aconvective core and a radiative envelope. Oscillation properties of the stellar model are also discussed in the context of observations.
Results.
Our tests show that specialized low-Mach techniques are necessary when simulating oscillations in stellar interiors. Basicproperties of internal gravity and pressure waves in our simulations are in good agreement with linear wave theory. As comparedto anelastic simulations of the same stellar model, we can follow internal gravity waves of much lower frequencies. The temporalfrequency spectra of velocity and temperature are flat and compatible with observed spectra of massive stars.
Conclusion.
The low-Mach compressible approach to hydrodynamical simulations of stellar oscillations is promising. Our simula-tions are less dissipative and require less luminosity boosting than comparable spectral simulations. The fully-compressible approachallows the coupling of gravity and pressure waves to be studied too.
Key words. hydrodynamics – methods: numerical – stars: interior – convection – waves
1. Introduction
The study of the excitation and propagation of waves withinstars has greatly helped to shape stellar structure and evolutiontheory over the last century. Today, this area of astronomy iscalled asteroseismology, and includes the study of oscillationsacross the Hertzsprung–Russell diagram. For example, pressuremodes (p-modes) provide important constraints on the envelopesof stars. Modes of consecutive radial order ( n ), and the sameangular degree ( (cid:96) ) have a characteristic frequency separationknown as the large frequency separation , which is sensitive tothe average density of a star (Aerts et al. 2010). This applica-tion of asteroseismology using p-modes has been extremely suc-cessful for low- and intermediate-mass stars (Chaplin & Miglio2013; Hekker & Christensen-Dalsgaard 2017; García & Ballot2019). Specifically, the measurement of envelope rotation usingrotationally-split pressure modes has facilitated the discoverythat stars with masses of about 2 M (cid:12) have approximately rigidinterior rotation profiles (see Kurtz et al. 2014; Saio et al. 2015;Van Reeth et al. 2016, 2018). Hence, current angular momen- tum theory needs significant improvement already on the mainsequence (MS) (Aerts et al. 2019).For later evolutionary stages, including subgiant, red giantand red clump stars, pulsation modes of a mixed character, i.e.pulsations that behave as gravity modes (g-modes) near the coreand as p-modes near the surface, have been detected in thousandsof stars (see Beck et al. 2011), and can be used to distinguishdi ff erent stages of nuclear burning (Bedding et al. 2011). Hence,understanding pressure modes is not only crucial for measuringinterior properties of main sequence stars, but also for post-MSstars (see e.g. Beck et al. 2012; Mosser et al. 2012).On the upper main sequence, stars with spectral types O andB ( M > (cid:12) ) observed in µ mag precision space photometryshow coherent opacity-driven p- or g-modes, as well as stochas-tic variability caused by internal gravity waves. This occurs inslowly pulsating B (SPB) stars with masses between 3 M (cid:12) and9 M (cid:12) (Pápics et al. 2017; Bowman et al. 2019a; Pedersen 2020),in β Cep stars with masses between 8 M (cid:12) and 25 M (cid:12) (Briquetet al. 2011; Burssens et al. 2019) and in young and evolved O-type dwarfs and blue supergiants with masses up to ∼
50 M (cid:12) (Buysschaert et al. 2015; Bowman et al. 2019b; Pedersen et al.
Article number, page 1 of 22 a r X i v : . [ a s t r o - ph . S R ] J un & A proofs: manuscript no. oscillations
Kepler / K2, and TESS photometry moti-vated the development and study of our simulation set up.The typical numerical approach to model the excitation ofwaves generated within stars is to solve the Navier–Stokes equa-tions in the anelastic approximation (e.g. Rogers et al. 2013;Alvan et al. 2014; Edelmann et al. 2019). This method allowslarge time steps while still being a mostly explicit method and istherefore computationally e ffi cient. While being suitable to sim-ulate internal gravity waves (IGWs) where gravity is the dom-inant restoring force, some important phenomena, such as theexcitation of p-modes, cannot be followed. Furthermore, com-mon numerical methods to solve the equations in the anelas-tic approximation require to introduce an artificial viscosity toachieve numerical stability. To balance the e ff ect of high viscos-ity, the stellar luminosity and the thermal di ff usivity have to beincreased by orders of magnitude. This leads to a damping of thewaves, especially in the low frequency regime.Some of these drawbacks are avoided or reduced by perform-ing compressible simulations of stellar interiors. Here, the fullNavier–Stokes equations are solved, most commonly in the finitevolume approach. This includes the physics of sound waves andallows the luminosity and the thermal di ff usivity to be kept muchcloser to stellar values. Viscosity is implicitly introduced by thenumerical scheme, and is lower than the viscosity typically usedby spectral codes (nevertheless still orders of magnitudes higherthan the astrophysical value). This is why we commonly speakof solving the Euler equations, i.e. the Navier–Stokes equationswithout an explicit viscosity term, in this context. On the otherhand, this kind of simulations comes with higher computationalcosts compared to their anelastic counterparts. That compress-ible simulations show excited IGWs and p-modes has alreadybeen reported in the past. For example, Meakin & Arnett (2006,2007) show a spectrum of the velocity for a simulation of carbonand oxygen burning in a 3D wedge geometry. They compare theform of a wave with predominately g-mode character and a cou-pled p- and g-mode to the predictions from linear theory and findgood agreement for both. Also Herwig et al. (2006) indicate theexcitation of p- and g-modes for the case of He-shell burning.They find the frequency of the di ff erent modes to be indepen-dent of resolution and boosting.These prior studies focus on the e ff ects of convective bound-ary mixing rather than the physics of waves. Thus, the compu-tational domains only contain small parts of the radiative zonebelow and above the convective shells in evolved stars. The fre-quency spectrum might therefore di ff er considerably comparedto a full star model and comparison to observations is di ffi cult.Also, the analysis of internal waves mainly consists of comput-ing the resulting spectra without further investigations.Therefore, to further assess the advantages of compress-ible simulations, we use our finite volume Seven-League Hy-dro (SLH) code (for a description see Sect. 2) to examine inmore detail the properties of excited waves. To ease the val-idation and comparison of our results, we repeat the simula-tion of a 3 M (cid:12) zero-age main-sequence (ZAMS) model by Edel-mann et al. (2019) (EM19 hereafter). For their simulation, EM19applied the anelastic approximation and a comparison betweenthese two approaches is therefore possible. We note that many ofthe diagnostics presented in the main part of this work originatefrom EM19.Because 3D simulations of an entire stellar model are costlyeven on today’s supercomputer facilities, we use 2D simulationsin this initial verification experiment. The 2D approach allows usto cover almost the entire stellar radius in our simulation domain, excluding only a small part in the core and the outermost layers,and to follow the evolution for an extended period of time. Thecomputational costs are low enough to run the simulations onthe local computer cluster of the
Heidelberg Institute for The-oretical Studies (HITS), Germany. Convection is known to beonly accurately described in three spatial dimensions (see e.g.Meakin & Arnett 2007 for a comparison of 2D and 3D oxygenburning). However, as pointed out by EM19, the results of their3D simulation are compatible with those of a 2D simulation byRogers et al. (2013) for a di ff erent, but similar 3 M (cid:12) model. Thisindicates that 2D simulations may still serve useful results forexcited wave properties despite their reduced dimensionality.This work is intended as a proof-of-concept: We present andvalidate in detail the results of simulating IGWs with the com-pressible hydrodynamics code SLH. The lack of the third di-mension may not allow a full quantitative comparison with ob-servations, but some characteristics of 2D simulations are stillexpected to match observations of the integrated variability atthe surface in a qualitative way. The main questions we aim toanswer are therefore: Can the excited waves be identified as pres-sure waves and IGWs and do they comply with the theoreticalexpectations? Additionally, we compare the excited wave spec-tra to previous 2D work by Rogers et al. (2013) and performqualitative comparisons to observations of frequency spectra ofmassive stars.The paper is organized as follows: In Section 2 we give abrief overview on the numerical methods that are used in theSLH code. In Section 3 we describe and apply a simple testsetup to benchmark the capability of SLH to treat IGWs. The2D results for the 3 M (cid:12) ZAMS model are discussed in Section 4.Section 5 summarizes the most important aspects and gives anoutlook for future simulations.
2. The SLH-Code
We use the SLH code for all simulations presented in this paper.It was developed initially by Miczek (2013) and solves the fullycompressible Euler equations in a finite volume framework. Itallows us to choose between an ideal gas and a more generalequation of state that includes contributions from radiation andelectron degeneracy (Timmes & Swesty 2000). The hydro solveris coupled to a nuclear reaction network (Edelmann 2014). Ra-diation is treated in the di ff usion limit, which is appropriate inthe optically thick regions of the star covered here. The imple-mented mapping procedure between the uniform Cartesian com-putational grid and a general curvilinear physical grid introducesflexibility regarding grid geometry. Our implementation of themapping is based on Kifonidis & Müller (2012) and Colella et al.(2011), examples for curvilinear grids can be found for examplein Calhoun et al. (2008).The SLH code was developed with a focus on flows at verylow Mach numbers. One complication in this regime is the needof specialized flux scheme as common approaches are domi-nated by artificial dissipation (e.g. see Miczek et al. 2015; Bar-sukow et al. 2017). For the simulations presented here, we usethe AUSM + -up solver (Liou 2006). It splits the flux function intoa pressure and an advective part and has improved low-Machcapabilities. An artificial di ff usive component for both parts isintroduced for stability. To avoid divergence at very low Machnumbers, the scaling of the di ff usion terms is limited by a cut-o ff Mach number
Mcut , which is a free parameter. For technicalreasons, a separate cut-o ff number Mcut pdi ff for the pressure dif-fusion term is used in SLH. In the simulations presented here,we use Mcut = − and Mcut pdi ff = .
1. For lower values of
Article number, page 2 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars
Mcut pdi ff , the convergence rate of the implicit solver decreasesconsiderably.Due to the large pressure gradient in stars, the AUSM + -upsolver quickly destroys the hydrostatic stratification. To resolvethis problem, SLH applies a variant of the Cargo–Leroux well-balancing technique. The basic idea is to remove the static back-ground stratification from the conserved variables before theyenter the numerical flux function. This considerably reducesthe “e ff ective” pressure gradient. The scheme was developed byCargo & Le Roux (1994) for one-dimensional setups. Edelmann(2014) describes how this approach can be extended to multi-dimensional setups.The flexible modular design of SLH facilitates the imple-mentation of new developments such that newly published fluxfunctions or reconstruction schemes can be easily implementedand tested. This makes SLH well suited to push hydrodynami-cal simulations towards low Mach numbers. For the 2D simu-lation of core hydrogen burning at the ZAMS presented in Sec-tion 4, mixing-length theory predicts convection at Mach num-bers around 10 − which calls for a numerical scheme optimizedfor slow flows. For a more detailed description of the code werefer the reader to Miczek (2013), Miczek et al. (2015) and Edel-mann (2014).At low Mach numbers, a common drawback of conventional explicit time stepping methods is that for stability the maximumpossible time step size δ t CFL uc is limited to the acoustic Courant–Friedrich–Lewy (CFL) criterion δ t CFL uc = CFL uc N dim min (cid:32) ∆ x | u | + c sound (cid:33) . (1)Here, N dim is the number of dimensions, ∆ x is the grid spac-ing in the di ff erent coordinate directions, u is the fluid velocityand c sound is the speed of sound. CFL uc is a dimensionless num-ber which needs to be smaller than unity. The factor 1 / N dim incombination with the minimum over all cells for the cell cross-ing time in all directions gives a lower limit for the time step.For low-Mach flows, u (cid:28) c sound and the time step size is domi-nated by the speed of sound. As a consequence, many small timesteps are needed to resolve the fluid flow. Although the compu-tational costs for one explicit time step are low, the large numberof necessary steps makes explicit time stepping ine ffi cient in thisregime.The great advantage of implicit time stepping is that there isno restriction for the step size required for stability. The mainconstraint arises from the question of how well the flow is to beresolved in time. This leads to the “advective” CFL u criterionwhich results from Eq. (1) by omitting the speed of sound: δ t CFL u = CFL u N dim min (cid:32) ∆ x | u | (cid:33) . (2)For illustration, setting CFL u = . Explicit first stage, Singly Diago-nally Implicit Runge-Kutta (ESDIRK) time steppers is imple-mented following the description of Hosea & Shampine (1996)and Kennedy & Carpenter (2001). For the simulations in this pa-per, the ESDIRK23 scheme is applied. It consists of three com-puting stages and is up to second order accurate in time. This results in two nonlinear systems of equations. SLH applies theNewton-Raphson method which finds their solution in an itera-tive way. Each iteration itself needs the solution of a linear sys-tem of equations. For these, SLH o ff ers a variety of di ff erentmethods (e.g. Krylov-Subspace schemes and multigrid solvers).At low Mach numbers, the increased computational costs pertime step are overcompensated by the benefit of a larger stepsize. Numerical tests with SLH imply that implicit time step-ping becomes more e ffi cient than explicit time stepping at Machnumbers smaller than about 0.1 to 0.01.Apart from accuracy requirements, the length of implicittime steps is also limited by the Newton-Raphson solver, whichmay converge slowly or even diverge if the time step becomestoo long. This limit depends on the problem solved and on de-tails of the numerical scheme and needs to be determined exper-imentally.The choice of the time step size for the 2D simulation pre-sented in Section 4 is discussed in Section 4.7 along with an e ffi -ciency comparison between implicit and explicit time stepping.
3. Testing internal gravity waves in SLH
In this section, we scrutinize the capability of SLH to correctlyreproduce the propagation of IGWs. The base setup is a 2DCartesian domain containing an isothermal ideal gas in a hy-drostatic stratification. A wave packet of small amplitude isevolved for several oscillation periods. The group velocity andthe change of the wave shape are then extracted from the simu-lation and compared to the theoretical prediction which followsfrom the linear Boussinesq approximation. With this simple butwell-defined test setup we verify whether SLH is able to repro-duce the prediction accurately enough before applying it to themore complex case of a realistic stellar profile in Section 4.
The theoretical basics of the benchmark test are presented in Ap-pendix A and follow Sutherland (2010). The actual experimentalSLH setup closely follows the idea of Miczek (2013) and is ex-tended up to Mach numbers of Ma = − . The basic parametersof the setup are listed in Table 1. Table 1.
Parameters of the Boussinesq IGW simulation. temperature T =
300 Kmean mol. weight µ = γ = / ρ = − gravity g = − cm s − e y resolution 288( x ) × y )domain [0 , (2 π/ | k x | )] × (cid:2) y , y (cid:3) boundary conditions x-direction: periodicy-direction: constant ghost cellstime stepping ESDIRK23time step size δ t =
120 2 π N reconstruction linear, second order in spaceAssuming the ideal gas law p = ρ R /µ T , where R is the uni-versal gas constant, the profiles for density, pressure, and poten- Article number, page 3 of 22 & A proofs: manuscript no. oscillations tial temperature in hydrostatic equilibrium are given by ρ hse = ρ exp (cid:32) − y H p (cid:33) , (3) p hse = ρ R µ T exp (cid:32) − y H p (cid:33) , (4) ϑ hse = T exp (cid:32) y H p γ − γ (cid:33) , (5)where the pressure scale height H p is defined as H − p = − ∂ ln p ∂y = gµ R T . (6)The Brunt–Väisälä frequency (BVF) according to Eq. (A.10) isspatially constant and reads N = (cid:115) g H p γ − γ . (7)We perturb this hydrostatic stratification with a monochromaticinternal gravity wave packet. The wavelength λ = π/ | k | = β H p of the packet is set in terms of the fraction β of the pressure scaleheight and is inclined by − ◦ with respect to the horizontal di-rection. We therefore have | k | = πβ H p , θ = − ◦ ◦ π (8)and k x = | k | cos θ, k y = | k | sin θ. (9)We set the vertical domain of the Cartesian box such that itcorresponds to 13 times the wavelength λy = − β H p , y = β H p (10)in order to provide enough space for the wave to move upwardsin y-direction. The horizontal extent is set to contain one hori-zontal wavelength λ x = π/ | k x | , which, in combination with pe-riodic boundary conditions, allows for the plane wave approach.At the top and bottom boundary we apply a layer of two cellswhich are filled with the hydrostatic initial condition but keptconstant in time (constant ghost cells). The amplitude for the ver-tical velocity component A v is modulated by a Gaussian functionaccording to A v ( x , t = = f Ma (cid:112) γ R T /µ (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) = c sound exp − (cid:32) yβ H p / (cid:33) . (11)The parameter f Ma therefore sets the peak Mach number inthe vertical velocity amplitude. Using the relations Eqs. (A.11)to (A.13) from the theory described in Appendix A one finds forthe initial conditions ϑ ( x , t = = ϑ hse + R (cid:40) − i ω d ϑ hse d y A v e i k · x (cid:41) , (12) u ( x , t = = R (cid:40) − k y k x A v e i k · x (cid:41) , (13) v ( x , t = = R (cid:110) A v e i k · x (cid:111) , (14) p ( x , t = = p hse + R (cid:40) − ρ ω k y k x A v e i k · x (cid:41) , (15) where R { . } denotes the real part of a complex expression.In Appendix A.1, the time evolution of an initial amplitudemodulation of the form of Eq. (11) is derived while consider-ing the IGW dispersion relation Eq. (A.9) up to second order in k . The result given by Eq. (A.28) shows that the initial profilebroadens in time and moves at a vertical velocity of c gy = | sin θ cos θ | N | k | (16)which corresponds to a vertical Mach number ofMa gy = | sin θ cos θ | (cid:112) γ − γ β π (17) ≈ . β. (18)In Section 3.2 we compare the predicted to the simulated evolu-tion of the velocity amplitude to assess the accuracy of our nu-merical schemes. In order to extract the velocity amplitude func-tion from the simulation, Miczek (2013) suggests the followingprocedure: – It is assumed that the horizontal velocity field can be decom-posed as u ( x , y ) = ˆ u ( y ) sin ( k x x + ϕ ( y )) . (19)This is fulfilled for the ansatz Eq. (A.6) and the ampli-tudes as defined above. Accordingly, also the evolved ini-tial data should at least approximately fulfill this decompo-sition. However, directly comparing this to the prediction isnot straight forward. – To simplify the interpretation, the square of Eq. (19) is inte-grated over the full horizontal width (cid:90) π/ k x u ( x , y ) d x = π k x ˆ u ( y ) . (20)The integral can be calculated numerically using the datafrom the simulation and therefore provides the possibility torecover the vertical profile ˆ u ( y ). – The normalized amplitude r j = (cid:113) k x π (cid:80) i u i , j ∆ x (cid:12)(cid:12)(cid:12)(cid:12) k y k x A v ( y = (cid:12)(cid:12)(cid:12)(cid:12) (21)then measures the change of the evolved data relative to themaximum of the initial data. In Eq. (21), i and j are the in-dices in horizontal and vertical direction, respectively; ∆ x refers to the size of the uniform grid spacing. – The velocity by which the peak of r j moves upwards is inter-preted as the group velocity in Eq. (16).As discussed in Section 4.6, one can define a nonlinearity pa-rameter, ε = u ω k x , (22)where u and k x denote the horizontal velocity and wave number,respectively. If ε (cid:38) ff ects to becomedominant. Inserting the corresponding expressions into Eq. (22)gives ε B = sin θ f Ma Ma gy (23)for the Boussinesq IGWs. Article number, page 4 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars
In this section we present simulations of the IGW setup de-scribed above for di ff erent parameter settings. For comparison,we show the results for the low-Mach solver AUSM + -up and theclassical Roe solver (Roe 1981).The possible parameter space is restricted by two conditions:(a) We require ε B (cid:28) β (cid:28) ε B = − and vary Ma gy . The val-ues for β and f Ma are then calculated accordingly. This way, weare able to assess the capabilities of both schemes for di ff erentMach numbers. The numerical settings for all of the simulationspresented in this section are listed in Table 1. In particular wehave chosen the time step size of the implicit time stepping suchthat the period corresponding to the BVF is resolved by 20 timesteps.We perform simulations for vertical group velocities ofMa gy = − , − , and 10 − . Whereas condition (b) is wellfulfilled for Ma gy = − , with Ma gy = − , which is closer tothe typical velocity we find the simulation presented in Section 4(see Fig. 4), we have β = × − and the stratification doesnot strictly follow the Boussinesq approximation anymore. Press(1981) shows that for the locally Boussinesq but globally anelas-tic equations the amplitude scales during the propagation with (cid:112) ρ /ρ where ρ is the density at the starting point (cf. Eq. (42)).We therefore multiply equation Eq. (A.28) by the correction fac-tor f ρ = (cid:118)(cid:116) ρ (cid:16) y − c gy t (cid:17) ρ ( y ) (24)to account for the amplification due to varying density.The results for the AUSM + -up and Roe solvers are visualizedin Figs. 1 to 3, respectively. The left columns show snapshotsof the horizontal velocity u at the end of the simulation. Theright columns compare at three points in time the amplitude andshape of the vertical velocity distribution as extracted from thesimulation using Eq. (21) with the approximate prediction givenby Eqs. (24) and (A.28).In Fig. 1, we have set Ma gy = − , which corresponds to avertical velocity amplitude of f Ma = − . For the AUSM + -upsolver the velocity field in the left column clearly shows the ef-fect dispersion: waves of longer vertical wavelengths move fasterand therefore appear at larger y values compared to smaller verti-cal wavelengths. The amplitude function broadens over time andis compatible with the prediction in terms of width and peak am-plitude. We attribute most of the small deviations from the pre-diction to our neglect of third order e ff ects in the dispersion re-lation, which are responsible for the amplitude function’s skew.In contrast, the Roe solver heavily damps the initial amplitudewithin the first few time steps and it becomes impossible to de-termine a unique peak in the remaining velocity field. This il-lustrates that the classical Roe solver fails in the very low Machregime whereas AUSM + -up still gives excellent results.We continue by increasing the vertical group velocity toMa gy = − . According to Eq. (17) this can only be achieved byincreasing β , which sets the fraction of the pressure scale heightcovered by one wavelength, to β = . × − . The results areshown in Fig. 2. The relative peak amplitudes and shapes forAUSM + -up do not change considerably compared to Fig. 1. Theresults of the Roe solver show distinguishable, but still strongly − . − . . . . . . y / H P u (cid:16) t = 6 · πN (cid:17) AUSM + -up Ma gy = 10 − , ε = 10 − f Ma ≈ − , β = 3 . × − t = 0 · πN t = 2 · πN t = 6 · πN .
000 0 .
003 0 . x/H P − . − . . . . . . y / H P . . . . . rel. amplitude of u x Roe
Fig. 1.
Results for the IGW test setup as described in Section 3 forthe AUSM + -up solver (upper row) and the classical Roe solver (lowerrow) . The parameters for the simulation are shown at the top of the plotand described in the main text. Left column:
Horizontal velocity u inthe 2D domain at the end of the simulation at t = π N . Blue color cor-responds to a positive value, and red color to a negative value of thevelocity. The scale is adjusted to the maximum amplitude for each run. Right column:
Amplitude extracted according to Eq. (21) at the begin-ning of the simulation and at two later points in time (solid lines). Theshaded areas correspond to the predicted shape of the amplitude mod-ulation function according to Eqs. (24) and (A.28). Dashed horizontallines mark the position of the peak amplitude for the prediction thatmoves at the group velocity according to Eq. (17). damped, peaks and their vertical group velocity considerablydisagrees with the theoretical prediction.The e ff ect of increasing the vertical group velocity further toMa gy = − is depicted in Fig. 3: Again, the shape and peak am-plitudes for AUSM + -up are compatible with the prediction. Forthis setup, β = × − and consequently the density varies no-ticeably in the simulation volume. This is reflected by the smallerdecrease of the expected peak amplitudes in Fig. 3 compared tothe amplitudes shown Figs. 1 and 2 for which the amplificationaccording to Eq. (24) is negligible. With the Roe solver, the re-sults are better as compared to those obtained at lower verticalgroup velocities and the broadening is closer to the prediction.However, significant damping of the velocity amplitude is stillevident.In summary, the results for the classical Roe solver, whichwe present here solely for comparison, clearly su ff er from highdamping and show that the wave packages move at the wrongspeed. Our tests therefore confirm the need for specialized low-Mach solvers in order to treat IGWs in the regime of velocitiesbelow Ma ∼ − . A variety of such solvers are readily availablein SLH. A promising example is the AUSM + -up solver for which Article number, page 5 of 22 & A proofs: manuscript no. oscillations − . − . . . . . . y / H P u (cid:16) t = 6 · πN (cid:17) AUSM + -up Ma gy = 10 − , ε = 10 − f Ma ≈ − , β = 3 . × − t = 0 · πN t = 2 · πN t = 6 · πN .
00 0 .
03 0 . x/H P − . − . . . . . . y / H P . . . . . rel. amplitude of u x Roe
Fig. 2.
See Fig. 1 for a description of the quantities shown. our tests demonstrate its capability to treat IGW in the low-Machregime.To assess the relevance of this finding we estimate the orderof magnitude of the group velocity we expect for the simulationof the real stellar setup presented in Section 4. We take the ab-solute value of Eq. (A.14) and rewrite it in terms of vertical andhorizontal components of the wave vector. For the polar coor-dinates used in the 2D simulation, the horizontal wave numbercorresponding to the angular degree (cid:96) is given by k h = √ (cid:96) ( (cid:96) + r (25)and the absolute value of the group velocity is | c g ( r ) | = r ω N √ (cid:96) ( (cid:96) + (cid:114) N ω − , (26)where r denotes the radial position within the model. In Fig. 4we show the group velocity for (cid:96) = − . These regions becomeeven more extended at higher (cid:96) -values. We conclude that low-Mach solvers are needed to correctly describe the wave field insimulations such as that presented in Section 4, especially in thelow-frequency regime dominated by internal gravity waves.
4. 2D simulation of a 3 M (cid:12) ZAMS star
The previous section demonstrates the capabilities of the SLHcode and the methods implemented therein to propagate IGWs inthe low-Mach regime whereas classical approaches fail. In this − . − . . . . . . y / H P u (cid:16) t = 6 · πN (cid:17) AUSM + -up Ma gy = 10 − , ε = 10 − f Ma ≈ − , β = 3 . × − t = 0 · πN t = 2 · πN t = 6 · πN . . . x/H P − . − . . . . . . y / H P . . . . . rel. amplitude of u x Roe
Fig. 3.
See Fig. 1 for a description of the quantities shown. f in µ Hz . . . . r / R ? ‘ = 4 − − − − − | v g | / c s o und Fig. 4.
Expected group velocity according to Eq. (26) expressed in termsof Mach number for values of the stellar model used in Section 4. Con-tour lines mark regions of di ff erent typical Mach numbers. There are noIGWs for frequencies above the BVF (white area). Instead, one expectsthe excitation of sound waves only which have Mach numbers of unity. section, the SLH code is applied to a real stellar setup, whichencompasses both the generation and propagation of IGWs andsound waves. For the simulation presented here we use the identical initial1D model as EM19. It describes a non-rotating 3 M (cid:12) star at theZAMS with a metallicity of Z = − and an outer radius ofR (cid:63) = . × cm. The model has been calculated with theopen-source stellar evolution code MESA (see Paxton et al. 2019for the latest report on code updates).The 1D data provided by the MESA model needs to bemapped onto the SLH grid while accurately fulfilling the equa- Article number, page 6 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars . . . . r/R ? − − − − ρ i n g c m − − − | N | π i n µ H z mappedMESA Fig. 5.
Density (blue) and absolute value of the BVF (orange). The dot-ted lines correspond to the profiles of the underlying 1D MESA model.The solid lines result from our mapping to the SLH grid. The verticaldashed line marks the radius below which the BVF has a negative signand one therefore expects convection to develop. The applied mappingmethod guarantees that the position of the dashed line is identical forthe mapped SLH model and the 1D MESA model. tion of hydrostatic equilibrium ∇ P = g ρ. (27)To do so, we re-integrate Eq. (27) while imposing the radial pro-file of one thermodynamic quantity from the 1D code (using the1D density profile for this would be a simple example). All otherthermodynamic quantities then follow from the equation of state.The particular choice of the imposed quantity depends on thespecific setup. For the case presented here, it is important to keepthe position of the convection zone as close as possible to the 1Dinput MESA model. Convective instability is characterized by anegative sign of the BVF. In regions without any or with only asmall gradient in composition (well fulfilled in ZAMS stars) it isessentially determined by the sign of ∇ − ∇ ad . Therefore, we fol-low the approach of Edelmann et al. (2017) to integrate Eq. (27)while enforcing( ∇ − ∇ ad ) SLH = ( ∇ − ∇ ad ) MESA . (28)This way the initial spatial extent of the convective zone on theSLH grid exactly matches the 1D input model. Consequently,other quantities might deviate from the 1D model. In Fig. 5 weexemplarily compare the profiles of the SLH model and the 1Dinput for density and BVF. Both quantities are reproduced rea-sonably well considering that we cover several orders of mag-nitudes.There is a deviation of a factor of ten in the convectionzone which we attribute to di ff erences in the technical details ofthe equation of state and the calculations of gradients betweenMESA and SLH. It is not expected that enforcing N = N will improve the result as this will only translate the di ff erencesinto other quantities. However, in the SLH simulation the valueof N in the convection zone will self-consistently adjust to avalue that corresponds to the equilibrium between energy in-put, e.g. due to nuclear burning, and the convective flux. Thesmall initial deviation is therefore not relevant. For the BVF, themean di ff erences in the radiation zone between the 1D MESAmodel and the mapped profiles is 0 .
56 % with a maximum valueof 1 . r = .
87 R (cid:63) and we therefore consider the appliedmethod to be su ffi ciently accurate.The underlying 1D MESA model is in thermal equilib-rium. The re-integration of Eq. (27) for hydrostatic equilibriumchanges the 1D stratification only slightly and we expect themapped SLH stratification in the radiative envelope to be suf-ficiently close to thermal equilibrium. This, however, is not truefor the convective core where we have to artificially boost thenuclear energy release (see Section 4.2). A stratification that is . . . . . . r/R ? κ ( r ) i n c m g − mappedMESA − − − τ d i ff ( δ r ) i n h Fig. 6. Blue lines: profiles of the opacity κ for the underlying 1D MESAmodel (dotted blue line) and the value applied for the mapped 2D SLHmodel (solid blue line). Orange lines: characteristic di ff usion time scaleaccording to Eq. (29) when assuming the radial spacing δ r of the SLHgrid as typical length scale. not in thermal equilibrium will re-adjust on the thermal-di ff usiontime scale τ di ff . It can be estimated as (e.g. Maeder 2009) τ di ff ( ∆ x di ff ) ∼ ( ∆ x di ff ) K , K = a c light T κ ρ C P , (29)where ∆ x di ff is a typical length scale. The thermal di ff usivity K introduces the radiation constant a = . × − erg cm − K − and the specific heat at constant pressure C P for the ideal gas.The opacity κ is a function of radius and determined by the phys-ical properties of the gas. However, for the mapped SLH model,we set κ SLH such that we achieve K MESA = K SLH . Therefore,the value of κ SLH is not fully consistent in a physical sense butit allows us to stay closer to the stellar values of thermal di ff u-sion. The blue lines in Fig. 6 show the profiles of κ SLH (solid)and κ MESA (dashed) and indicate that deviations are reasonablysmall.Following Eq. (29), we estimate the thermal-di ff usion timescale taking the radial grid spacing δ r to be the typical lengthscale, see Fig. 6. Due to the energy boosting in the convectionzone, deviations from thermal equilibrium are expected to belargest there. From Fig. 6 it can be seen that the di ff usion timescale is of the order of 10 h and thus two orders of magnitudelarger than what is covered by our 2D simulation (about 10 h,see Section 4.2). Furthermore, as the grid spacing δ r is alreadythe smallest possible scale, our estimate is a lower limit on thetime scale. In the outer parts, where the time scale becomes com-parable to the simulation time, our model is expected to be closeenough to thermal equilibrium.As EM19 do not need to modify the 1D input MESA model,these plots also serve as a direct comparison between the mappedmodel on the SLH grid and their initial data.In Table 2 we list the parameters for the SLH simulation. Weuse 2D polar geometry which corresponds to an infinite cylin-der. SLH currently does not support the geometry of a 3D spherewith azimuthal / longitudinal symmetry. Because of the singular-ity of polar coordinates at r =
0, we cannot include the wholecore. The minimum radius of the domain is mainly determinedby the decreasing cell sizes in horizontal direction which af-fects the possible time step size according to the CFL criterion(see Section 4.7). We have chosen r min = .
007 R (cid:63) which stillallows reasonably large steps. The upper boundary was set to r max = .
91 R (cid:63) which is close to the value of EM19 who set theupper boundary to 0 . (cid:63) . We apply solid-wall boundary con-ditions at the inner and outer boundaries of the computationaldomain. They enforce a vanishing velocity perpendicular to theboundary interface. This prohibits mass flux and sound wavesfrom leaving the domain. Periodic boundaries are chosen in the Article number, page 7 of 22 & A proofs: manuscript no. oscillations
Table 2.
Parameters of the 2D simulation. R (cid:63) = . × cm denotesthe total radius of the underlying 1D model, L MESA the stellar luminosityand K MESA the stellar thermal di ff usivity as given by MESA. boosting L = L MESA th. di ff usivity K = K MESA
EoS ideal gas + radiation pressuregeometry polar coordinatesradial domain 0 .
007 R (cid:63) to 0 .
912 R (cid:63) resolution 960( r ) × ϕ )boundary conditions r -direction: solid-wall ϕ -direction: periodictime stepping ESDIRK23time step size δ t = + -upazimuthal direction, which is appropriate since we cover the fullazimuthal range of 2 π . The number of 960 radial cells ensuresthat the smallest pressure scale height (close to the outer bound-ary) is still resolved by 16 grid cells. The number of horizontalgrid cells is set such that the cell width at the top of the convec-tion zone δ w = δ ϕ r top roughly matches the height δ r of the cells,where δ ϕ = π/ N ϕ and δ r = ( r max − r min ) / N r denote the angularand radial resolution, respectively. For the stellar luminosity as given in the 1D MESA input model,mixing-length theory (MLT, e.g. Kippenhahn et al. 2012) pre-dicts Mach numbers for the convective core of around Ma ∼ − . Simulations in this regime are numerically very challeng-ing and the SLH code has several specialized low-Mach approxi-mate Riemann solvers implemented. However, we have recentlynoticed that for Mach numbers considerably below 10 − the con-vective flow seems not to be driven by heating but rather bynumerical instabilities. This problem is subject of ongoing in-vestigations but there is no adequate solution available yet. As aworkaround, it was decided to artificially boost the energy gener-ation in order to increase the convective velocity. From MLT onefinds that the convective velocity v conv scales with the luminosity L as v conv ∝ L / (30)which has also been confirmed by numerous numerical studies(e.g. figure 7 of Cristini et al. 2019 or figure 15 of Andrassy et al.2019). We boost the stellar luminosity by a factor of 10 whichcorresponds to a tenfold increase in velocity compared to theMLT prediction. We note that this boosting is still a factor 10 smaller compared to the boosting in the simulations of EM19and Rogers et al. (2013).In Fig. 7 we compare the MLT prediction to the results ofour 2D simulation. The gray shaded area at small radii marksthe region of the core that is convectively unstable according tothe Schwarzschild criterion. In the input 1D model, an additionalsmall convection zone near the surface of the star is present butour 2D model already ends at a smaller radius. We find our sim-ulations in good agreement with this “scaled” MLT prediction.However, outside of the convective core a velocity field ofconsiderable amplitude has developed. These velocities are at-tributed to the excitation and propagation of IGWs and are of . . . . . . r/R ? − − − − − M a c h N u m b e r v rms , v ∗ conv , MLT v conv , MLT
Fig. 7.
Predicted and simulated Mach numbers of the 2D model. The or-ange lines correspond to the MLT prediction: The dashed profile showsthe original stellar values, whereas the solid line illustrates the veloc-ities scaled by 10, according to Eq. (30). The blue line shows the 2Dresults averaged over roughly two convective turn over times, startingat t =
500 h. The profile ends at r ≈ . R (cid:63) as our domain does not con-tain the entire star model. The gray shaded areas mark the convectivecore and the small region of surface convection in the 1D model. main interest for the work presented here. In stark contrast tothat, no velocities are assumed in the input 1D model. This il-lustrates the shortcomings of 1D stellar evolution, where thesedynamical phenomena have to be parametrized (see e.g. Aertset al. 2019).The speed of sound in our model ranges from 7 × cm s − within the convection zone to 1 × cm s − at the top of thecomputational domain. Accordingly, the Mach number shown inFig. 7 corresponds to typical velocities of 2 × cm s − (con-vection zone) to 1 . × cm s − (near surface). For the meanconvective turnover time in the convection zone τ conv we find τ conv = ∆ r cz v rms ≈
40 h (31)where ∆ r rc = . × cm is the depth of the convection zoneand v rms the root-mean-square of the absolute velocity. The spa-tial mean has been taken over the convection zone and the tem-poral mean includes the entire simulation except for the initialtransient phase, where convection has not yet developed (see alsoFig. 10). Hence, the 700 h of fully developed convection that wefollow in our simulation cover roughly 17 convective turnovertimes.The velocity field of the 2D simulation is further illustratedin Fig. 8. It visualizes radial velocity (left panel) and tempera-ture fluctuations from the horizontal average (right panel). Bothquantities are scaled by their horizontal root-mean-square as theamplitudes vary considerably between the outer and inner part ofthe model. The relative amplitudes have maximum values of fourand are approximately equal at all radii for temperature and ve-locity. In both panels of Fig. 8, the convective core is filled witha few rather coherent structures which correspond to convectiveeddies that induce wave-like motions. These waves form spiralpaths from the point of excitation at the boundary of the convec-tive core as they travel towards the surface. The spatial structuresare smaller than those observed by EM19 (e.g. see their figure 8).This is explained by the much smaller e ff ective viscosity of oursimulation and further illustrated in Fig. 9 which shows an SLHsimulation of a similar domain and resolution but with explicitviscosity and thermal di ff usivity comparable to those used byEM19. It is clearly visible how the enhanced di ff usive e ff ectseven out small scale structures and only large patterns remain.Figs. 7 and 8 illustrate that the fundamental process of gen-erating internal waves, i.e. a convective core with plumes thatexcite waves in the radiative envelope, is present in our 2D sim- Article number, page 8 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars − . − . − .
25 0 .
00 0 .
25 0 .
50 0 . x/R ? − . − . − . . . . . y / R ? v r − . − . − .
25 0 .
00 0 .
25 0 .
50 0 . x/R ? T fluc − − q ( r , ϕ ) / q ( r ) r m s Fig. 8.
Radial velocity v r and temperature fluctuations T fluc . As both quantities have larger amplitudes at the outer parts of the model (see furtherSection 4.5) they have been scaled with the corresponding horizontal root-mean-square value at each radius r to ease the visualization. Both panelsalso show a magnification of the core region. The shown snapshot is taken at t ∼
580 h. There is also a movie available. − . − . − . . . . . x/R ? − . − . − . . . . . y / R ? v r − . − . − . . . . . x/R ? T fluc − − q ( r , ϕ ) / q ( r ) r m s Fig. 9.
Same quantities as in Fig. 8 but here for a 2D simulation with a viscosity and thermal di ff usivity similar to the values used by EM19. Notethat the domain for this simulation is comparable to Fig. 8 but not identical. The radial resolution is slightly higher for the viscous simulation anddue to the high viscosity the energy input is boosted by 10 in order to get convection starting. ulation. In order to further validate our results in the subsequentsections, we closely follow the methods as presented in EM19. The fundamental properties of stellar oscillations can be probedby investigating their temporal frequency spectrum. For this, weperform a temporal Fourier transformation (FT) of the entire ve-locity field in our 2D simulation.The transformation is done for both the radial and horizon-tal velocity. In order to select waves corresponding to a specificangular degree (cid:96) , we apply a filter prior to the temporal FT. For the velocity of each stored snapshot, a spatial FT in the angle ϕ is taken while using all available data points on the grid. Subse-quently, all amplitudes are set to zero, except for the one corre-sponding to the (cid:96) value we want to filter for. This manipulatedspectrum is brought back to real space via an inverse FT. The re-sulting time sequence of the grid cells of each radial ray is thenmultiplied by the Hanning window to reduce leakage e ff ects andserves as input for the temporal FT. To reduce the backgroundnoise such that individual modes appear more clearly, the tempo-ral spectra are taken for 100 evenly distributed radial rays. Thesquares of the amplitudes of the resulting Fourier coe ffi cientsare then averaged. In principle, the average could be done for all Article number, page 9 of 22 & A proofs: manuscript no. oscillations available angles but we experienced no improvement for a largernumber of rays. Keeping the numbers of rays low is desirableregarding memory requirements.The coe ffi cients of the temporal FT are normalized in thesame way as in EM19 (see their equation 12 and 13), i.e. the am-plitudes are divided by the number of the input data points. Thisresults in coe ffi cients that are independent of the number of binsand that have the same units as the input quantity. Furthermore,for this normalization, the amplitude of a peak across one singlefrequency bin corresponds to the actual amplitude of that wavein the time domain which eases the direct interpretation of thespectrum. For our spectra, the frequency bins have a width of δ f = / (700 h) = . µ Hz. Spectral lines in our spectra typically(except for a few broad p-modes) span two to four frequency binsand the absolute amplitude is therefore slightly underestimatedwhen directly read o ff the spectral plots. Finally, the amplitudesare multiplied by a factor of 2 to account for the change in am-plitudes in the frequency spectrum due to the application of the Hanning window (e.g. Sect. 9.3.9 of Brandt 2011). t in h − − − M a c h N u m b e r Ma max Ma rms Fig. 10.
Maximum (blue) and root-mean-square (orange) Mach numberin the CZ as a function of time. The gray shaded area marks the timeframe that has been used to extract the spectra that are presented in thispaper.
We only use the time interval spanning 700 h from 400 hto 1100 h for the spectral analysis (see the gray shaded area inFig. 10) to avoid the initial development of convection in thecore and of the wave field in the envelope.In Fig. 11 we show the result for the radial velocity, whichis the dominant component for p-modes. The upper panel showsthe spectrum of the unfiltered velocity, the panels below showthe spectra of filtered velocities for some exemplary values of (cid:96) .This figure is the compressible counterpart of figures 24 and 27in EM19. In Fig. 12, the same quantities are shown for the hor-izontal velocity which is the dominant component for g-modes.In these plots and in the further course of the paper, a hat denotesquantities that have been obtained using a temporal FT.The white line in Fig. 11 indicates the profile of the BVF atthe start of the simulation, the black dashed line corresponds tothe profile averaged over the last 100 hours of the simulation. Amagnified vision of the top of the domain is given for (cid:96) = ff ect is likely relatedto the boundary conditions and grid resolution, the exact causeis not fully understood at this point and still subject of ongoingwork. However, the deviation is only within a small part of thewhole model and we do not expect it to have a significant impacton the general results nor the frequency spectrum of the wavesin the envelope below R = . (cid:63) . To resolve also the detailsof the outermost parts more accurately, a more elaborate outerboundary condition is needed as well as a higher resolution toresolve the small scale heights accurately enough. IGWs can only exist for frequencies below the BVF (e.g. seeSec. 3.4.2 in Aerts et al. 2010). This property is reflected in thefirst panel of Fig. 11 by the fact that the whole possible IGWregime (as defined by f < N / π ) shows a significant amplitudewhich then suddenly drops for f > N / π . This is in qualitativeagreement with figure 27 of EM19 and a clear indicator of theexcitation of IGWs in the 2D SLH simulation. An increasingnumber of radial nodes (or, equivalently, a decrease of the radialwavelength) for lower frequencies is another general property ofIGW that is confirmed by the spectra in Fig. 11.These features have already been seen and described inEM19. However, our simulation also shows additional signals.For (cid:96) = (cid:96) > S (cid:96) = (cid:96) ( (cid:96) + c r , (32)(e.g. Sect. 3.3.2 of Aerts et al. 2010). They are derived for spher-ical geometry using spherical harmonics and might not strictlyhold in our 2D geometry. However, we still use them here as anestimate to characterize the general behavior of internal waves.A p-mode of angular degree (cid:96) and frequency f p can only exist ifthe frequency fulfills the conditions f p > N π and f p > S (cid:96) π . (33)The corresponding conditions for g-modes are f g < N π and f g < S (cid:96) π , (34)(see Sec. 3.4 of Aerts et al. 2010 for further details). In regionswhere the frequency does not fulfill these relations, p- and g-modes are evanescent. This is most easily seen for p-modes inthe spectrum for (cid:96) = (cid:96) =
20 in Fig. 11.The occurrence of p-modes is one of the key di ff erences be-tween fully compressible codes and the anelastic approaches. Inthe latter, p-modes cannot occur because sound waves are notincluded in the equations. To bring the p-modes of our 2D simu-lation into context with observations, we compare them to thoseof the β Cep stars presented in Aerts & De Cat (2003). They findfrequencies at low (cid:96) typically around f ≈ − ≈ µ Hz. Forthe particular star ω Sco observations indicate f =
15 d − ≈ µ Hz at (cid:96) =
9. Of course, the eigenfrequencies of a real stardepend on its stellar parameters and excitation mechanism, and β Cep stars have masses typically in the range of 8 M (cid:12) to 20 M (cid:12) .Thus, the BVFs which set the minimum p-mode frequencies (seeEq. (33)) are smaller compared to our 3 M (cid:12) model. We find am-plitudes associated with p-modes starting at around 300 µ Hz.Although our model is not directly comparable to the observed β Cep stars, this indicates that the waves, which are excited in oursimulations by the convective core, are compatible with thoseobserved in real stars of similar mass. We note, however, thatthis cannot be seen as a proof for the correctness of the model orthe underlying excitation mechanism, especially since p -modesare not typically observed in main-sequence 3 M (cid:12) stars (Aertset al. 2010).
Article number, page 10 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars . . . . . r / R ? all ‘N/ π S S S S S S . . . . . r / R ? ‘ = 00 . . . . . r / R ? ‘ = 10 . . . . . r / R ? ‘ = 20 . . . . . r / R ? ‘ = 30 . . . . . r / R ? ‘ = 40 . . . . . r / R ? ‘ = 50 200 400 600 800 1000 f in µ Hz . . . . . r / R ? ‘ = 20 10 − − − ˆ v r i n c m s − Fig. 11.
Frequency spectrum of the radial velocity for a time span of 700 h. The normalization is done in a way that the amplitude of a narrow lineof one frequency bin width corresponds to the velocity amplitude of the corresponding wave in the time domain. The data is stored in intervalsof 480 s. This allows us to capture frequencies up to f max = µ Hz with a resolution of δ f = . µ Hz. In order to reduce the background noise,we show the average of the spectrum of 100 individual radial rays. The doublets for the modes with f (cid:38) µ Hz are due to aliasing, which weverified with a simulation with a very short cadence of outputs. The three colored dots in the uppermost panel mark the radii for which the lineprofiles are shown in Fig. 16. The white and black lines correspond to the BVF at start and end of the simulation, respectively. Colored dashedlines correspond to the Lamb frequencies of di ff erent (cid:96) values according to Eq. (32). For (cid:96) =
0, the uppermost part of the model is magnified (graybox) in order to illustrate the change in the BVF at the end of the simulation. For the magnification, the color scale was adapted to increase thevisibility of the lines. We show the spectra for the horizontal velocity in Fig. 12. Article number, page 11 of 22 & A proofs: manuscript no. oscillations . . . . . r / R ? all ‘N/ π S S S S S S . . . . . r / R ? ‘ = 00 . . . . . r / R ? ‘ = 10 . . . . . r / R ? ‘ = 20 . . . . . r / R ? ‘ = 30 . . . . . r / R ? ‘ = 40 . . . . . r / R ? ‘ = 50 200 400 600 800 1000 f in µ Hz . . . . . r / R ? ‘ = 20 10 − − − ˆ v h i n c m s − Fig. 12.
Frequency spectrum for the horizontal velocity, see Fig. 11 for details. The p-modes are less prominent as their dominant velocitycomponent is radial.
For radii deep inside the convection zone ( r (cid:46) .
14 R (cid:63) ) dis-tinct modes are less pronounced in Figs. 11 and 12 (g-modesgenerally cannot exist there and in our case the Lamb frequenciesprohibit also p-modes) but instead amplitudes in a wide range offrequencies appear, a consequence of the turbulent convection inthe core. This is most easily seen for larger (cid:96) values in the spec-trum of the horizontal velocity. On the other hand, also at the largest radii, a band of high amplitudes is visible in the upper-most panel which extends over all frequencies. This is due to thedevelopment of a shear flow towards the end of the simulation.As it is located very close to the outer boundary, we cannot de-termine beyond doubt whether the shear flow is caused by theboundary or by deposition of the IGW’s angular momentum due
Article number, page 12 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars to nonlinear e ff ects. This needs to be examined in more detail infuture simulations.Another way of illustrating the emerging spectra of waves isgiven in Fig. 13 (this is similar to figure 22 in Herwig et al. 2006or figure 2 in Alvan et al. 2015). Here, the spectrum at a fixedradius r = . R (cid:63) is shown for the first 100 (cid:96) -values. The hori-zontal blue line indicates the BVF. One clearly observes distinctmodes for small (cid:96) and decreasing mode spacing for increasing (cid:96) .The amplitudes drop at the BVF as expected for gravity waves.Alvan et al. (2015) estimate a cut-o ff frequency ω c below whichwaves have lost a considerable amount of the kinetic energy dueto radiative di ff usion (see their equation 26). The correspondingprofile is shown as dashed white line. They argue that a travelingwave without su ffi cient energy cannot reach the turning point,travel back and interfere with another progressive wave to forma standing mode. Their corresponding cut-o ff is rather steep andreaches high frequencies, indicating a wide region where theyexpect progressive waves. This is similar to figure 5 of Rogerset al. (2013). In contrast, for the simulation presented here theprofile of ω c is much flatter and stays at much smaller frequen-cies. It is clear that the excitation of gravity waves from core con-vection produces an entire spectrum of waves spanning a broadrange in frequency.In Fig. 14 we plot the kinetic energies v ( (cid:96) ) = (cid:88) f ˆ v (cid:96), f , v ( ω ) = (cid:88) (cid:96) ˆ v (cid:96),ω (35)at the top of the convection zone. Here, ˆ v (cid:96),ω is the velocity atangular degree (cid:96) and frequency ω which results from the spa-tial and temporal FT. The sum over all degrees (cid:96) is terminatedat (cid:96) = v ( (cid:96) ), shown in our left panel.Also the position of the transition between the two regimes iscomparable. However, our profiles for v ( ω ) (right panel) aremuch steeper (Rogers et al. 2013 find ω − . and ω − . ). In the 2Dsimulation presented here the transition between the two regimesis at much higher frequencies. So far we have presented clear indications of the excitation ofIGWs and pressure waves. To further confirm this, we show inthe following section that the corresponding dispersion relationsare fulfilled for both g- and p-modes.
From theory, one finds a dispersion relation for IGWs of the form(e.g. Aerts et al. 2010, Sec. 3.1.4 or Sutherland 2010, Sec. 3.3.3) ω N = k h | k | with | k | = (cid:113) k h + k r , (36)where ω = π f is the angular frequency of the gravity waveand k h , k r are the horizontal and radial wave numbers. We notethat this expression is derived under the assumption of a spatiallyconstant value for the BVF. From Fig. 5 it is clear that this is notthe case for the stellar model presented here. Therefore, Eq. (36)is expected to hold only for wavelengths that are short comparedto the scale height of the BVF.To estimate how close our simulation follows Eq. (36), weapply the method of EM19 which is briefly summarized in thefollowing: The angular frequency ω and BVF are readily availablequantities. Because of the spatial FT filtering for specific (cid:96) , weknow that k h = (cid:96) ( (cid:96) + / r . The only quantity to determine inEq. (36) is therefore the vertical wave number k r = π/λ r orequivalently the radial wavelength λ r . To determine the wave-length, the positions of peaks in the amplitude for the FT of theradial velocity are identified for all available frequencies. Thedistance between adjacent peaks is interpreted as one half of thewavelength at the radial position at the midpoint between thetwo peaks. We determine the wavelengths starting just above theconvection zone until the upper boundary of the model for allfrequencies and interpolate in radius afterwards. This gives foreach frequency f the wavelength λ r , f ( r ). From that, we obtain k r and are finally able to evaluate Eq. (36) in the entire radius-frequency plane for a specific order (cid:96) .To visualize the results, the ratio of the left and right side ofEq. (36) is plotted in Fig. 15 exemplarily for (cid:96) =
2, 4, 10 and 20.Regions in which the dispersion relation is fulfilled, i.e. the ratiois unity, appear white. This reproduces figure 28 of EM19. A ma-genta line denotes the frequency below which one wavelength isresolved by less than 10 grid points: For a given radial grid spac-ing δ r , the corresponding wave number is k r , n = = π/ (10 δ r )such that this frequency is defined as f n = = N π k h (cid:113) k h + k r , n = . (37)The yellow line in Fig. 15 marks the frequency below whichone assumes the waves to be dissipated by thermal di ff usion.The estimate is based on the equality of the di ff usion lengthand the wavelength at the critical frequency as given by equa-tion
27 of EM19. In comparison to EM19, the e ff ect of dis-sipation is greatly reduced in the SLH simulation. In our case,we are mainly limited by the radial resolution when going tohigher (cid:96) -values. We see agreement with the dispersion relation oflinear internal gravity waves everywhere where such agreementcan be expected, i.e. for modes well resolved in space but withradial wavelengths short enough that N as well as the radiuscan be considered approximately constant over a single wave-length. This is the case for low-frequency IGW and gets worsefor increasing radial wavelengths, i.e. increasing frequencies. Incontrast to EM19, our results are mainly a ff ected by resolutione ff ects rather than damping, yet broad spectra of standing modesare clearly excited in both numerical setups. For a better visibil-ity, the low frequency regions are magnified in the left columnof Fig. 15. At low frequencies (i.e. below 5 µ Hz), it is di ffi cult todetermine if the dispersion relation is fulfilled because the qual-ity of the FT deteriorates as the number of wave periods that fitinto the time window of our analysis decreases. However, we areable to reach much smaller frequencies than the anelastic simu-lations and provide reliable results for wave frequencies above5 µ Hz.Another way of testing the dispersion relation is to measurethe inclination of the IGW crests with respect to the radial direc-tion. As the fluid velocity in an internal wave is parallel to thewave crest, we can obtain the inclination by measuring the ratio v h /v r in our simulation. Because the wave vector k is perpendic-ular to the crests we have (cid:32) v r v h (cid:33) (cid:107) (cid:32) − k h k r (cid:33) , where (cid:32) − k h k r (cid:33) · k = . (38) Their equation is missing a factor (cid:112) N /ω −
1. It does not changethe result qualitatively but is included here for consistency.Article number, page 13 of 22 & A proofs: manuscript no. oscillations angular degree ‘ f i n µ H z r = 0 . R ? − − ˆ v h i n c m s − Fig. 13.
Color map of amplitudes at di ff erent frequencies as a function of the angular degree (cid:96) at a fixed radius of R = . (cid:63) . The blue dashedline corresponds to the BVF at this radius. The white dashed lines indicates the cut-o ff frequency ω c below which waves are expected to lose toomuch kinetic energy due to damping to form standing modes. ‘ ˆ v i n c m s − ‘ − . ‘ − . ω in µ rad s − − ˆ v i n c m s − ω − . ω − . r = 0 . R ? Fig. 14.
Kinetic energy as a function of angular degree (cid:96) (left panel)and angular frequency ω (right panel) at the top of the convection zone.Dashed lines correspond to power-law fits. This figure is similar to thesecond and third column in figure 6 of Rogers et al. (2013). It follows that (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v h v r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) k r k h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (39)and with the aid of Eq. (36) we find the expression (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v h v r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = √ N − ω ω , (40)in which the left-hand side can be obtained from the simula-tion and compared to the theoretical prediction on the right-handside.In the top two panels of Fig. 16 we show a line plot rep-resentation for the spectra of v r and v h at three di ff erent radialpositions. They are marked by colored dots in the upper panelof Fig. 11. For radii well above the convection zone, the am-plitudes rapidly decrease for f > N / π . This is expected for asignal mostly made up of IGWs which cannot propagate in thisfrequency range. Above frequencies of N / π we can see clearlyisolated peaks corresponding to p-modes of various angular de-grees. The p-mode amplitudes are at least one order of magni-tude smaller than the IGW signal across a broad range of fre-quencies. For the radial velocity, the spectrum is almost flat witha subtle decrease towards higher frequencies. The spectrum forthe horizontal velocity shows a decrease already within the IGWregime. Similar to the findings of EM19, these integrated spec-tra composed of all angular degrees do not show distinct peaks. However, our results do not have the steep drop in amplitude atvery low frequencies as seen in their figure 23, which can beattributed to the lower viscosity and thermal di ff usivity in oursimulation.In the lowest panel of Fig. 16 we show the ratio of the ve-locity components as depicted in the two panels above (semi-transparent lines). The solid lines represent the prediction givenby Eq. (40). There is an excellent agreement for all three radialpositions for frequencies below the BVF. The line for 0 . R (cid:63) is just beyond the boundary of the convective core and thereforedoes not follow the steep drop in the ratio as convective motionsimpact the data. Only for the smallest frequencies the simulationdeviates as a result of damping e ff ects and a lack of independentdata points for the temporal FT. This result is another strongindication that IGWs are correctly represented in the 2D simula-tion.Furthermore, we find the results of our simulations to becomparable to observations of late-B SPB stars. De Cat & Aerts(2002) report the ratio v h /v r for several SPB stars which are typ-ical g-mode pulsators (this corresponds to the K -value in theirfigure 17). The value ranges from 10 to 100. This observed rangeis similar to the values shown in the lower panel of Fig. 16 forfrequencies to the left of the dip, i.e. below the BVF. For p-modepulsators De Cat (2002) (see also Aerts & De Cat 2003) reports K -values for β Cep stars which are typically in the range of 0 . (cid:96) ) are suppressedin this quantity as they cancel out on average whereas large-scalecontributions are pronounced (e.g. Aerts et al. 2010, Sect. 6). Asa simple estimate of this e ff ect we compute the average temper-ature over half of the circumference of our 2D model prior to thetemporal FT. The result is shown as an orange line in Fig. 17 (forcomparison, we show as a blue line the spectra which have been Article number, page 14 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars ‘ = 20 . . . . r / R ? µ Hz ‘ = 40 . . . . r / R ? µ Hz ‘ = 100 . . . . r / R ? µ Hz f in µ Hz ‘ = 20 f max . . . . r / R ? µ Hz − ( k h / k ) / ( ω / N ) Fig. 15.
Ratio of the left and right side of the dispersion relation Eq. (36).The method to extract k h / k from the simulation is described in the text. Aratio of unity corresponds to a perfect match and is reflected by white regions. A red color indicates a too large value of k h / k and correspondinglythe ratio is too small in blue regions. The white line is the BVF. For frequencies below the yellow line, IGWs are expected to be damped bythermal di ff usion. Below the magenta line, the radial wavelength is resolved by less than 10 grid cells. The left column zooms into the lowfrequency region of the corresponding full plane on the right. The frequency given in the white boxes denotes the value of f max which sets the scaleof the corresponding x-axis. averaged similar to the velocity spectra, i.e. only after the tem-poral FT). At high frequencies, pressure modes appear clearlyas individual peaks. This is, however, not the case in the low fre-quency regime. We note that our spectrum is like those presentedin Bowman et al. (2019a) (e.g. see their figure 3 for the observ-ables of the O star HD 46150 and the figures in their appendixfor additional OB stars). The actual situation for observationsis more complex than the diagnostic value of our simulationscan be, e.g. we do not account for limb darkening, while thisphenomenon comes into play to determining the net flux varia-tion. The limb darkening has more e ff ect for flux observationsbut far less for the net velocity variations in the line-of-sight.However, our simple experiment illustrates that the net flux vari-ations do not easily reveal any distinct peaks corresponding tolow- (cid:96) modes in the low-frequency regime despite the expectedcancellation of the numerous high- (cid:96) modes. For the pressure waves, we perform a similar analysis as de-scribed in the previous section except that we now check for theusual dispersion relation of pressure waves ω c sound = | k | . (41)The result is shown in Fig. 18 and the colors have the samemeaning as in Fig. 15. As expected, the dispersion relation is notfulfilled for f < N / π . At frequencies corresponding to a stand-ing mode we get a good agreement with predictions from theory.Furthermore, Fig. 18 indicates that pressure modes are excited in the mixed-mode frequency regime, i.e. between about 200 µ Hzand 300 µ Hz. This suggests that a coupling of p- and g-modes ispossible. The deviation between standing modes is probably dueto the same reasons as discussed in Section 4.4.1. This is furtherillustrated in Fig. B.1 of the Appendix where we show the ampli-tudes at two specific frequencies. When the frequency matchesa standing mode the wavelength can be detected easily whereasthis is not possible for frequencies in-between. These results to-gether with the fact that the modes for a particular angular degreeare equidistant in frequency in the regime of high-order p-modesgive us confidence that the observed modes are indeed p-modes.
From linear theory, one expects the amplification of IGW to-wards the surface due to density stratification. As shown in Rat-nasingam et al. (2019) for spherical geometry, the prediction forIGW amplification is given by (in the notation of EM19) v h ∝ (cid:18) r r (cid:19) / (cid:114) ρ ρ N − ω N − ω / exp( − τ/ , (42)where r is the starting point of the wave and ρ , N the corre-sponding density and BVF. The damping factor τ is given by τ = (cid:90) rr d r κ ( (cid:96) ( (cid:96) + / N r ω (cid:114) − ω N . (43) We note that there is an error in the corresponding equation 23 ofEM19. The expression holds for the horizontal velocity component v h instead of the vertical velocity v r as stated.Article number, page 15 of 22 & A proofs: manuscript no. oscillations − − ˆ v h i n c m s − f o r a ll ‘ fit function: a · f b b = − . , r = 0 . R ? b = − . , r = 0 . R ? b = − . , r = 0 . R ? − − ˆ v r i n c m s − f o r a ll ‘ b = − . , r = 0 . R ? b = − . , r = 0 . R ? b = − . , r = 0 . R ? f in µ Hz − ˆ v h / ˆ v r − f in d − Fig. 16.
Line plot representations of the spectra for horizontal and verti-cal velocity at three di ff erent radii 0 .
14 R (cid:63) , 0 . (cid:63) and 0 . (cid:63) are shownin the upper two panels . The dashed lines correspond to the power lawfit as denoted in the labels. The lowest panel shows the ratio of the ve-locity components at each radius (transparent line) and compares themto the prediction according to Eq. (40) (solid lines). f in µ Hz − − ˆ T i n K r = 0.80 R ? all ‘ all ‘ mean ϕ ∈ [0 , π ] Fit a · f b , b = − . BVF − f in d − Fig. 17.
Spectra for temperature at r = .
80 R (cid:63) . The orange line showsthe spectrum of the temperature fluctuations which have been averagedover half of the circumference of our 2D model. The blue line corre-sponds to the spectra averaged over 100 radial rays after the FT as it isdone for the velocity spectra. The dashed blue vertical line denotes theposition of the BVF. The gray line corresponds to a power-law fit for thefrequencies ranging from 10 µ Hz to 200 µ Hz resulting in an exponent of b ≈ − According to Eq. (42), waves are damped by the e ff ects ofthermal dissipation and geometry, but amplified by decreasingdensity. This ignores the damping e ff ect of viscosity. For con-sistency, Eq. (42) needs to be multiplied by a correction fac-tor ( r / r ) − / to account for our 2D polar geometry of an in-finite cylinder, which slightly increases the amplification (Rat- . . . r / R ? ‘ = 00 . . . r / R ? ‘ = 20 200 400 600 800 1000 f in µ Hz . . . r / R ? ‘ = 4 10 − | ~ k | ω / c s o und Fig. 18.
Same quantities as in Fig. 15 but this time for the dispersionrelation of sound waves Eq. (41). ˆ v h i n c m s − f‘ =1 = 6 . µ Hz f‘ =4 = 13 . µ Hz f‘ =1 = 76 . µ Hz MESA extrapolations . . . . . . r/R ? ˆ v r i n c m s − Fig. 19.
Velocity amplitude for di ff erent frequencies and angular de-grees (cid:96) as function of radius. The amplitudes also include the contri-bution of the two adjacent frequency bins to account for the fact thatpeaks typically show a width of two to four bins in our simulation. Theexpected amplification towards the surface according to Eq. (42) is rep-resented by dashed lines; the 1D MESA profile data serve as input vari-ables whereas the starting points are taken as the simulated amplitudesat the first noticeable local maximum of the respective frequency andangular degree. nasingam et al. 2020, submitted). For the radial velocity, Eq. (42)needs to be scaled according to Eq. (40).In Fig. 19 we compare the result of Eq. (42) for the 3 M (cid:12) MESA model to the 2D simulation, similar to figure 26 ofEM19. For both radial and horizontal velocity, the predictionfrom linear theory and the simulation are in good agreement forshort radial wavelengths and low frequencies ( f (cid:96) = = . µ Hz, f (cid:96) = = . µ Hz) and we assume that the MESA extrapola-tions towards the surface provide a reasonable estimate of sur-face velocities. At the highest frequency and longest wavelength( f (cid:96) = = . µ Hz) prediction and simulation clearly di ff er. How-ever, this is expected as the prediction assumes that wavelengthsare short compared to all relevant scale heights in the stratifica-tion which is clearly not the case for the highest frequency shownin Fig. 19.Our extrapolation towards larger radii is a very simplifiedapproach and we do not account for the complex physics at thesurface, such as the existence of a sub-surface convection zone. Article number, page 16 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars . . . . . . r/R ? − − − − − − n o n li n e a r i t y p a r a m e t e r ε f‘ =1 = 6 . µ Hz f‘ =4 = 13 . µ Hz f‘ =1 = 76 . µ Hz MESA extrapolations
Fig. 20.
Nonlinearity parameter ε according to Eq. (44) for the samefrequencies and angular degrees as in Fig. 19. Therefore, these numbers should only be seen as an approxima-tion to the order of magnitude of the velocity at the surface dueto amplified stellar oscillations. The drop in the simulation dataat the very top of the computational domain is an artifact of thenumerical solid-wall boundary condition whereas the drop in theprediction is an e ff ect of stronger radiative damping near the sur-face.Also for the extrapolated velocities at the stellar surface, theorder of magnitude for the ratio of the horizontal velocity overthe radial one shown in Fig. 19 is compatible with what is foundin time-series spectroscopy of g-mode pulsators (De Cat & Aerts2002). If all waves were to remain in the linear regime, a treatment asin Section 4.5 would be a su ffi cient description of the physics in-volved. Yet it is expected that the amplification towards the sur-face causes nonlinearities to become relevant at least in certainranges of frequency and wavenumber. These nonlinearities leadto a redistribution of energy between di ff erent wavenumbers andfrequencies. Additionally, wave breaking can cause transport ofangular momentum from the core, where the waves are excited,to the envelope.In Ratnasingam et al. (2019) the e ff ect of di ff erent inputspectra on the expected nonlinearity of IGWs is examined. Thenonlinearity can be estimated by ε = v h ω k h (44)(Press 1981; Barker & Ogilvie 2010). If ε (cid:38)
1, nonlinear e ff ectsare dominant. However, as demonstrated in Ratnasingam et al.(2020, submitted), already rather small values of ε ( ≈ − ) cancause noticeable energy transfer between di ff erent frequenciesand wavenumbers.In Fig. 20 we show the result of Eq. (44) for the amplitudesof the horizontal velocity extracted from the simulation and theextrapolation to the surface as illustrated in Fig. 19. The apparentincrease in ε with decreasing frequency is caused by the strongerconvective wave excitation at lower frequencies. At even lowerfrequencies, which are not shown here, wave damping becomesdominant and ε decreases further. For the extrapolated ampli-tudes we find a maximum value for the nonlinearity parameterof ε = . f (cid:96) = = . µ Hz.From this we conclude that nonlinear e ff ects can be expectedat the surface in the frequency regime around 10 µ Hz. This isfurther illustrated in Fig. 21, which shows the value for ε at r = . (cid:63) as a function of frequency f and angular degree (cid:96) .A narrow range around 10 µ Hz is apparent where the nonlinear- f in µ Hz ‘ r = 0 . R ? − − − n o n li n e a r i t y p a r a m e t e r ε Fig. 21.
Nonlinearity parameter at r = . (cid:63) for all available frequen-cies up to an angular degree of (cid:96) = . . . . r/R ? δ t i n s δt CFLuc=0 . δt CFLuc=50 δt CFLu=0 . δt sim Fig. 22.
Time step size as a function of radius for di ff erent time stepcriteria. Because of the polar geometry of the grid, the innermost cellsbecome narrower and thus require shorter time steps. The purple dashedhorizontal line denotes the time step size of δ t = ity is highest. We note that this looks very similar to SpectrumK and
Spectrum LD in figure 5 of Ratnasingam et al. (2019) forconvective velocity boosted by a factor of three with respect tothe stellar models. The strongest nonlinearity in Fig. 21 is seen inthe frequency range in which gravity modes have been detectedin a sample of about 30 SPB stars in the
Kepler data. For all thesepulsators, nonlinear behavior has been deduced from the
Kepler light curves and a low-frequency IGW power excess has beendetected after the removal of the coherent g-modes (Pedersen2020).We emphasize that the predictive power for the wave am-plification and nonlinearity from this SLH simulation is limited.The boosted energy generation in the core and the accompany-ing increased velocities also impact the amplitudes of the ex-cited waves. Furthermore, convection is known to be faster in 2Dsimulations and we apply a simple extrapolation from the upperboundary of our 2D model towards the surface. Nevertheless,we believe that the results presented here at least give an idea ofwhat to expect at the stellar surface where we find indication forthe existence of nonlinear e ff ects. After the detailed evaluation of the properties of the excitedwaves in the previous sections, the paper is concluded with adiscussion on the time step sizes and the e ffi ciency of the SLHcode.In Fig. 22, the radial profiles of the time step size δ t CFL uc = . ( r )for explicit time stepping (see Eq. (1)) and the maximum im-plicit time step size δ t CFL u = . ( r ) (see Eq. (2)) are shown for val-ues in the middle of the time span of the 2D simulation. Thetime step sizes decrease towards the core due to the shrinking Article number, page 17 of 22 & A proofs: manuscript no. oscillations azimuthal width of the cells in polar geometry. From the radialprofiles, the actual global time step is then given by the mini-mum over the whole domain. We find min[ δ t CFL uc = . ( r )] = .
03 sand min[ δ t CFL u = . ( r )] = . δ t ≤ δ t CFL uc = , whereasthey get considerably damped but still propagate at a speed closeto the theoretical prediction for δ t ≈ δ t CFL uc = . For δ t ≥ δ t CFL uc = sound waves cannot be followed at all anymore.For the 2D simulation presented here, these findings indicatethat choosing a time step size of δ t CFL u = . δ t ≈ δ t CFL uc = at the lower boundary might lead to a strongdamping of the sound waves. Although the appropriate time stepsize clearly depends on the details of the setup and the wave-lengths of the considered waves, we follow Miczek (2013) inchoosing the time step size in this initial study. As will be shownat the end of this section, this might give a value for the time stepsize that is too conservative.For the particular simulation presented here, the time stepsize is chosen to be δ t sim = ff er from damping as for these cells the timestep is larger than δ t CFL uc = . This, however, is only a tiny fractionof the whole domain and we do not expect it to have a signifi-cant impact on the results. Our choice is a compromise betweene ffi ciency and accuracy as will be demonstrated in the following.To quantify the gain in e ffi ciency when using implicit timestepping, we compare the wall clock time needed to cover a timespan of 80 min (2 convective turnover times) when using eitherimplicit or explicit time stepping. For this test, both simulationsare restarted at 500 h and all output routines of SLH are dis-abled in order to minimize possible external e ff ects like a slowfile system. We find that the explicit run requires a wall clocktime of 37 .
72 min to perform 15 900 steps on 360 Intel Skylakecores. The implicit run finishes after 6 . ffi ciency of the SLH code (implicitand explicit) to the pseudo-spectral anelastic SPIN code (EM19).For the comparison, the test runs presented in this section areused for the SLH timings. For the same physical problem, a shorttest simulation on 300 Ivy Bridge cores for a grid of 1500( r ) × ϑ ) × ϕ ) serves as input for the SPIN results.In Table 3, the timings for the two codes to perform one timestep per cell and core are listed. From this measure, the explicitSLH mode is roughly a factor of five slower than the SPIN code,potentially reflecting the fact that SLH has not undergone anysubstantial optimization e ff ort. The implicit SLH mode is 220times slower than SPIN with a time step size that is only eighttimes larger. These numbers illustrates that, while our compress- ible simulations capture the full set of physics described by thefull Euler equations, they come with considerably higher com-putational costs.We note that the settings of the SLH implicit mode maynot be optimal yet. As stated by Miczek (2013), fine tuning thelinear-solver for the specific physical problem to be solved andapplying a multigrid solver may lead to improved performance.This is subject of ongoing work and the necessary algorithms arereadily available in SLH. Table 3.
Timings for the SLH code ( implicit and explicit ) and thepseudo-spectra anelastic code SPIN to perform one time step per celland core.
Code time / cell / step (cid:2) w − µ s (cid:3) SLH (impl) 3 . × SLH (expl) 7 . . JUWELS su-percomputer in Jülich, Germany. The domain is discretized in280( r ) × ϑ ) × ϕ ) cells while performing time steps of δ t = × core-h will beneeded to cover 700 h of simulation time. By scaling this numberwith the number of grid cells, a simulation with 480 × × . × core-h. The reference sim-ulation applies the same time step size as the finer resolved 2Dsimulation; therefore, a possible change in the step size is notincluded in the scaling. However, depending on the desired ac-curacy, a much larger time step can be chosen and the computa-tional costs decrease correspondingly.This demand of computational resources is a realistic sce-nario for applications at HPC facilities. A 3D simulation at thesame resolution as the 2D simulation requires 44 × core-hwhich is at the upper limit of common computing time propos-als. However, these numbers show that SLH is e ffi cient enoughto perform simulations of stellar oscillations also in 3D for suf-ficiently long time spans.As demonstrated, the implicit time stepping improves the ef-ficiency of the SLH code. At the same time it is important to con-firm that the accuracy at which sound waves are evolved is su ffi -cient. To this end, we restart the simulation from the 2D implicitrun at 500 h and evolve it using an explicit, three-stage Runge-Kutta (RK3) scheme (Shu & Osher 1988) which is third-orderaccurate in time. The time step is set to δ t CFL uc = . , our standardchoice for explicit time stepping. The simulation is evolved for100 h of physical time which corresponds to roughly 80 soundcrossing times t cross . We define the sound crossing time to be thetime a sound wave needs to travel from the innermost point ofthe computational domain to the upper boundary, t cross = (cid:90) r r c sound ( r ) d r = .
26 h . (45)We expect the time span of 80 t cross to be su ffi cient to reveal pos-sible deviations in the spectra. Article number, page 18 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars − f = 320 . µ Hz , ‘ = 4 RK3 δt ≈ .
03 s
ESDIRK23 δt = 8 s10 − − f = 1005 . µ Hz , ‘ = 1 − f = 280 . µ Hz , ‘ = 1 . . . . . r/R ? f = 41 . µ Hz , ‘ = 4 . . . . . . . . . . . . ˆ v r i n c m s − Fig. 23.
Amplitudes for the radial velocity spectra at four di ff erent fre-quencies for explicit (blue) and implicit (dashed orange) time stepping.The first two frequencies correspond to p-modes whereas the lower twofrequencies match g-modes. For radii above the convection zone, theamplitudes completely overlap each other. In Fig. 23, the spectra of the implicit and explicit run arecompared at di ff erent frequencies corresponding to g- and p-modes. All spectra are based on the same time frame.Almost no di ff erences between the amplitudes from implicitand explicit time stepping are visible. This result confirms thatthe propagation of sound waves is treated correctly also in theimplicit run. Furthermore, our implicit time step size of 8 s ap-pears to be a rather conservative choice and larger steps sizesmight be possible. This could reduce the computational costsconsiderably and will be investigated in more detail in futuresimulations.
5. Conclusion
The main goal of this paper has been to verify the capabilitiesof the time-implicit, compressible SLH code to correctly treatIGWs and p-modes in stars with a convective core. To this end,two test cases have been considered.The first test is a simple 2D setup where an IGW packet ac-cording to the Boussinesq approximation is evolved in time.We first simulated the propagation of an IGW packet in aweakly-stratified 2D atmosphere. The initial IGW was set upwith di ff erent expected group velocities. Velocity distributionswere extracted from simulations performed using the low-MachAUSM + -up solver and the classical Roe solver at a few points intime. The ability of the two solvers to follow IGW was assessedby comparing the numerical solutions with approximate analyticsolutions. The Roe solver revealed strong damping and di ff er-ent broadening of the initial wave packet and the packet’s groupvelocity was incorrect. These e ff ects were most pronounced atMach numbers Ma gy = − . The AUSM + -up scheme showedno significant damping and propagated the packet at a speedclose to the prediction. This indicates that a specialized solveris necessary when treating IGWs at low Mach numbers.Our second test case involving low-Mach-number flows wascore hydrogen burning in a 3 M (cid:12) ZAMS star. We used the sameinitial 1D model as EM19 to make the simulations comparable. A 2D simulation of this setup was performed using the AUSM + -up solver. The properties of IGW and p-modes were studied fol-lowing similar methods as in EM19. It was shown that the spec-tra extracted from the 2D SLH simulations reflect the basic prop-erties of internal waves.A broad spectrum of IGWs is observed for the integratedspectrum whereas individual standing modes can be identified inspectra for single angular degrees (cid:96) . Modes below the BVF havean increasing radial order with decreasing frequency, a funda-mental property of IGWs. Also the dispersion relation extractedfrom the simulation and the ratios of vertical to horizontal ve-locities match the theoretical prediction. Furthermore, we findthe velocity ratios to be compatible with observational diagnos-tics from time-series space photometry and high-resolution spec-troscopy of SPB stars (De Cat & Aerts 2002).For standing modes above the BVF, the radial order in-creases with increasing frequency and the dispersion relationmatches the one of p-modes. This kind of waves cannot be seenin the anelastic approximation as sound waves are removed fromthe underlying equations.Recently, Lecoanet et al. (2019) argued that if the observedvariability of stellar surfaces was due to the excitation of IGWfrom core convection, one would expect to observe distinct peaksin the spectrum which correspond to low (cid:96) values (see their fig-ure 2). We do not see such features in our simulation (e.g. seeFig. 17). Our explanation for the broad and near-flat profile ofIGWs is the same as in EM19: The small spacing between reso-nant frequencies of di ff erent radial order and di ff erent (cid:96) “hides”individual peaks. Moreover, stellar rotation, which is not in-cluded in the simulations presented here, would cause spectralline splitting and a further increase in the number of lines in thespectrum.The amplification of the waves towards the surface agreeswith the expectation given in Ratnasingam et al. (2019). In SLH,the treatment of IGWs at very low frequencies is mainly limitedby radial resolution. In contrast, anelastic simulations are lim-ited by radiative damping and viscous e ff ects already at largerfrequencies. Irrespective of the numerical setup, this work andthat of EM19 demonstrate the importance of the excitation andpropagation of IGWs as a diagnostic tool for the interior physicsof stars burning hydrogen in a convective core.The simulation of the 3 M (cid:12) model presented here is intendedas a proof of concept and aids in the comparison of the sim-ulations of Rogers et al. (2013) and Edelmann et al. (2019).The chosen 2D geometry reduces computational costs and al-lows for parameter exploration. A validation of 2D results basedon selected 3D models is planned for future work. From thecost of preliminary low-resolution 3D simulations we estimateda need of 44 × core-h to simulate a grid with a size of960( r ) × ϑ ) × ϕ ) for 700 h physical time. Using onlyhalf of the number of cells in each dimension reduces the esti-mated cost to 5 . × core-h. These estimates are based on areference run with the same time step size as the higher resolved2D simulation presented here. Thus, a change in the time stepsize is not considered in the scaling. Our tests indicate that theimplicit time step size could be increased while still resolvingsound waves accurately enough. This could considerably reducethe computational costs. Finally, SLH has already proven an ex-cellent scaling up to a large number of cores (e.g. Edelmann &Röpke 2016) such that this kind 3D simulations are feasible onmodern HPC facilities.After having tested the methods on the 3 M (cid:12) model, we willextend our study to higher stellar masses and later evolution-ary stages for which there are more observational data. Further- Article number, page 19 of 22 & A proofs: manuscript no. oscillations more, we aim to use the velocity and temporal information fromour hydrodynamics data to extract synthetic observables by aver-aging appropriately over the di ff erent wavenumber components.Studying the dependence of wave amplitudes on di ff erent lumi-nosity boosting will help us to estimate the amplitudes by ex-trapolating towards stellar values. Acknowledgements.
LH, FKR, and RA acknowledge support by the KlausTschira Foundation. PVFE was supported by STFC grant ST / L005549 / References
Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology(Springer)Aerts, C. & De Cat, P. 2003, Space Sci. Rev., 105, 453Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, 57, 35Alvan, L., Brun, A. S., & Mathis, S. 2014, A&A, 565, A42Alvan, L., Strugarek, A., Brun, A. S., Mathis, S., & Garcia, R. A. 2015, A&A,581, A112Andrassy, R., Herwig, F., Woodward, P., & Ritter, C. 2019, MNRAS, 2556Barker, A. J. & Ogilvie, G. I. 2010, MNRAS, 404, 1849Barsukow, W., Edelmann, P. V. F., Klingenberg, C., Miczek, F., & Röpke, F. K.2017, Journal of Scientific Computing, 72, 623Beck, P. G., Bedding, T. R., Mosser, B., et al. 2011, Science, 332, 205Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55Bedding, T. R., Mosser, B., Huber, D., et al. 2011, Nature, 471, 608Bowman, D. M., Aerts, C., Johnston, C., et al. 2019a, A&A, 621, A135Bowman, D. M., Burssens, S., Pedersen, M. G., et al. 2019b, Nature Astronomy,3, 760Bowman, D. M., Burssens, S., Símon Díaz, S., et al. 2020, A&A, acceptedBrandt, A. 2011, Noise and Vibration Analysis: Signal Analysis and Experimen-tal Procedures, EBL-Schweitzer (Wiley)Briquet, M., Aerts, C., Baglin, A., et al. 2011, A&A, 527, A112Burssens, S., Bowman, D. M., Aerts, C., et al. 2019, MNRAS, 489, 1304Buysschaert, B., Aerts, C., Bloemen, S., et al. 2015, MNRAS, 453, 89Calhoun, D. A., Helzel, C., & Leveque, R. J. 2008, SIAM Review, 50, 723Cargo, P. & Le Roux, A. 1994, Comptes rendus de l’Académie des sciences.Série 1, Mathématique, 318, 73Chaplin, W. J. & Miglio, A. 2013, ARA&A, 51, 353Colella, P., Dorr, M. R., Hittinger, J. A., & Martin, D. F. 2011, Journal of Com-putational Physics, 230, 2952Cristini, A., Hirschi, R., Meakin, C., et al. 2019, MNRAS, 484, 4645De Cat, P. 2002, Astronomical Society of the Pacific Conference Series, Vol. 259,An Observational Overview of Pulsations in β Cep Stars and Slowly Pulsat-ing B Stars (invited paper), ed. C. Aerts, T. R. Bedding, & J. Christensen-Dalsgaard, 196De Cat, P. & Aerts, C. 2002, A&A, 393, 965Edelmann, P. V. F. 2014, Dissertation, Technische Universität MünchenEdelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4Edelmann, P. V. F. & Röpke, F. K. 2016, in JUQUEEN Extreme Scaling Work-shop 2016, ed. D. Brömmel, W. Frings, & B. J. N. Wylie, JSC Internal ReportNo. FZJ-JSC-IB-2016-01, 63–67Edelmann, P. V. F., Röpke, F. K., Hirschi, R., Georgy, C., & Jones, S. 2017,A&A, 604, A25García, R. A. & Ballot, J. 2019, Living Reviews in Solar Physics, 16, 4Hekker, S. & Christensen-Dalsgaard, J. 2017, A&A Rev., 25, 1Herwig, F., Freytag, B., Hueckstaedt, R. M., & Timmes, F. X. 2006, ApJ, 642,1057Hosea, M. & Shampine, L. 1996, Applied Numerical Mathematics, 20, 21 ,method of Lines for Time-Dependent ProblemsKennedy, C. A. & Carpenter, M. H. 2001, Additive Runge-Kutta Schemesfor Convection-Di ff usion-Reaction Equations, Tech. rep., NASA TechnicalMemorandum Kifonidis, K. & Müller, E. 2012, A&A, 544, A47Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution(Berlin Heidelberg: Springer-Verlag)Kurtz, D. W., Saio, H., Takata, M., et al. 2014, MNRAS, 444, 102Lecoanet, D., Cantiello, M., Quataert, E., et al. 2019, ApJ, 886, L15Liou, M.-S. 2006, Journal of Computational Physics, 214, 137Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars, Astron-omy and Astrophysics Library (Springer Berlin Heidelberg)Meakin, C. A. & Arnett, D. 2006, ApJ, 637, L53Meakin, C. A. & Arnett, D. 2007, ApJ, 667, 448Miczek, F. 2013, Dissertation, Technische Universität MünchenMiczek, F., Röpke, F. K., & Edelmann, P. V. F. 2015, A&A, 576, A50Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10Pápics, P. I., Tkachenko, A., Van Reeth, T., et al. 2017, A&A, 598, A74Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10Pedersen, M. G. 2020, Phd thesis, KU Leuven in Belgium, in pressPedersen, M. G., Chowdhury, S., Johnston, C., et al. 2019, ApJ, 872, L9Press, W. H. 1981, ApJ, 245, 286Ratnasingam, R. P., Edelmann, P. V. F., & Rogers, T. M. 2019, MNRAS, 482,5500Ratnasingam, R. P., Edelmann, P. V. F., & Rogers, T. M. 2020, submitted toMNRASRoe, P. L. 1981, Journal of Computational Physics, 43, 357Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772,21Saio, H., Kurtz, D. W., Takata, M., et al. 2015, MNRAS, 447, 3264Shu, C.-W. & Osher, S. 1988, Journal of Computational Physics, 77, 439Sutherland, B. R. 2010, Internal gravity waves (Cambridge university press)Timmes, F. X. & Swesty, F. D. 2000, ApJS, 126, 501Van Reeth, T., Mombarg, J. S. G., Mathis, S., et al. 2018, A&A, 618, A24Van Reeth, T., Tkachenko, A., & Aerts, C. 2016, A&A, 593, A120 Article number, page 20 of 22orst et al.: Fully compressible simulations of waves and core convection in main-sequence stars
Appendix A: Linear theory in Boussinesqapproximation
In order to get analytical expressions describing the behavior ofIGWs, the fully compressible Euler equations need to be lin-earized. This can be done in the Boussinesq approximation. Itis based on the assumption that pressure and density vary onlylittle over in the volume considered. Furthermore, one imposes atime-independent hydrostatic background state and only followsthe temporal evolution of deviations from the background state.Additionally, a divergence-free velocity field is assumed. Thisapproach removes the physics of sound waves. For a Boussinesqgas it is convenient to introduce potential temperature ϑ as ϑ = T (cid:32) pp (cid:33) − ( γ − /γ , (A.1)where p is the pressure at a specific reference height in the at-mosphere with ϑ , ρ (see Sutherland 2010 for a detailed intro-duction). The two-dimensional equations of motions can then bewritten asD ˜ ϑ D t + v d ϑ hse d y = , (A.2)D u D t + ρ ∂ ˜ p ∂ x = , (A.3)D v D t + ρ ∂ ˜ p ∂y = − gϑ ˜ ϑ, (A.4) ∂ u ∂ x + ∂v∂y = , (A.5)where quantities with a tilde denote fluctuations from the hy-drostatic background state, e.g. ˜ p = p − p hse . The letters u , v refer to the horizontal and vertical components of the velocity u = ( v, u ) T . Further, the notation above makes use of the materialderivative D q / D t = ∂ q /∂ t + u ∇ q .A solution to this set of equations can be found using theansatz of a 2D plane wave ˜ ϑ u v ˜ p = A ϑ A u A v A p exp [ i ( k · x ) − i ω t ] , (A.6)which introduces the wave vector k = ( k x , k y ) T , the angular ve-locity ω , and the complex amplitudes A i . The absolute values ofthese amplitudes are assumed to be small, such that terms withproducts of two or more amplitudes can be neglected. This es-sentially removes the advection term in the material derivative.Inserting the ansatz into Eqs. (A.2) to (A.5) results in a homoge-neous system of linear equations of the form − i ω d ϑ hse d y − i ω ik x ρ gϑ − i ω ik y ρ ik x ik y · A ϑ A u A v A p = M · A = . (A.7)Non-trivial solutions exist only if det( M ) = ω = − d ϑ hse d y gϑ k x | k | = − d ϑ hse d y gϑ cos ( θ ) , (A.8)where we have used k · e x = k x = | k | cos θ with e x being the unitvector in horizontal direction and θ the angle between the wave vector k and the horizontal direction. The dispersion relation isusually written in the form ω = N cos( θ ) , (A.9)where N = (cid:115) − d ϑ hse d y gϑ hse ≈ (cid:115) − d ϑ hse d y gϑ (A.10)is the BVF in the Boussinesq approximation (e.g. see Suther-land 2010, Sec. 3.2). This result shows that the angular frequencydoes not depend on the absolute value of the wave vector k andthat IGWs do not propagate isotropically. The maximum fre-quency is ω = N for θ = ◦ and no purely vertical waves existas ω = θ = ◦ .For the specific solution of Eq. (A.7), we set the amplitudeof the vertical velocity as free parameter and express the otheramplitudes accordingly: A ϑ = − i ω d ϑ hse d y A v , (A.11) A u = − k y k x A v , (A.12) A p = − ρ ω k y k x A v . (A.13)As Eq. (A.7) is a linear system, any superposition of solutionsremains a solution to the system. The group velocity of such awave packet is then given by c g = ∇ k ω = N k x cos θ sin θ (cid:32) sin θ − cos θ (cid:33) , (A.14)if one uses sin θ = k y / | k | . Appendix A.1: Evolution of a wave packet
The following section is based on the methods described inSutherland (2010, Sect. 1.15) which we generalize to 2D andapply to the specific setup presented in Section 3.In linear theory, a quasi-monochromatic wave packet η ( x , t )is usually described as η ( x , t ) = A ( x , t ) exp [ i ( k · x − ω ( k ) t )] , (A.15)where the amplitude modulation function A ( x , t ) changes muchslower in space and time than the exponential function inEq. (A.15) and can be represented via a FT as η ( x , t ) = (cid:90) ∞−∞ ˆ η ( k ) exp [ i ( k · x − ω ( k ) t )] d k . (A.16)Applying the inverse FT to Eq. (A.16) and using Eq. (A.15), weobtainˆ η ( k ) = π ) (cid:90) ∞−∞ η ( x ,
0) exp [ − i k · x ] d x (A.17)and furthermore A ( x , t ) = η ( x , t ) exp [ − i ( k · x − ω ( k ))] = (cid:90) ∞−∞ ˆ η ( k ) exp [ i ∆ k · x ] exp [ − i ∆ ω t ] d k , (A.18) Article number, page 21 of 22 & A proofs: manuscript no. oscillations where ∆ k = k − k and ∆ ω = ω ( k ) − ω ( k ). For t =
0, Eq. (A.18)illustrates that an initial amplitude modulation A ( x ,
0) intro-duces waves with k (cid:44) k . Because the wave packet is quasi-monochromatic with a typical wavenumber of k , the amplitudeˆ η ( k ) must decrease quickly for k (cid:44) k . However, in the presenceof a nontrivial dispersion relation ω = ω ( k ), any superpositionof waves with di ff erent wavenumbers will lead to a broadeningof the initial wave packet over time. This is will be explicitlyshown for the setup presented in Section 3.Here, the initial vertical velocity modulation is a Gaussian asgiven by Eq. (11): η ( x , = A exp (cid:34) − y σ (cid:35)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) A ( x , exp [ i k · x ] (A.19)with A = f Ma (cid:112) γ R T /µ ; σ = β H p / . (A.20)Evaluating the FT in Eq. (A.16) leads toˆ η ( k ) = δ ( k x − k x , ) √ πσ A exp (cid:34) − σ ∆ k y (cid:35) . (A.21)where δ denotes the Dirac delta function. This specific form of η ( k ) illustrates that a narrower Gaussian modulation in real spaceleads to a broader distribution in wavenumber space and conse-quently to a larger dispersion.From Eqs. (A.18) and (A.21) the time evolution of A is de-termined. To evaluate the FT in Eq. (A.18), we expand the dis-persion relation in a Taylor series up to second order ω ( k ) = ω ( k ) + ∇ k ω | k ∆ k + (cid:16) ∂ k x , k x ω | k ∆ k x + ∂ k y , k y ω | k ∆ k y (cid:17) + ∂ k x , k y ω | k ∆ k x ∆ k y + O (cid:16) k (cid:17) (A.22)where we have introduced the abbreviations ∂ q = ∂∂ q , ∂ p , q = ∂ ∂ p ∂ q for convenience. Inserting Eqs. (A.21) and (A.22) intoEq. (A.18) gives A ( x , t ) = √ πσ A (cid:90) ∞−∞ exp (cid:104) i ∆ k y y (cid:105) · exp (cid:34) − σ ∆ k y − i (cid:16) ∇ k ω | k + ∂ k y , k y ω | k ∆ k y (cid:17) t (cid:35) d k y . (A.23)This can be simplified to A ( X , t ) = √ πσ A (cid:90) ∞−∞ exp (cid:34) − ˜ σ ∆ k y (cid:35) exp (cid:104) i ∆ k y Y (cid:105) d k y , (A.24)by transforming into a moving frame via X = ( X , Y ) T = x − ∇ k ω | k t (A.25)and setting˜ σ = (cid:113) σ + i ∂ k y , k y ω | k t . (A.26)The integral in Eq. (A.24) is again the FT of a Gaussian and leadsto the final expression A ( X , t ) = A σ ˜ σ exp (cid:34) − Y σ (cid:35) (A.27) where the absolute value of A is given by |A ( X , t ) | = A σ (cid:32) σ + (cid:18) ∂ k y , k y ω | k t (cid:19) (cid:33) / · exp − Y σ + (cid:32) ∂ k y, k y ω | k t σ (cid:33) . (A.28)Eqs. (A.27) and (A.28) describe the broadening and decrease inamplitude of an initial Gaussian profile over time when consid-ering the dispersion relation up to second order in k in a systemof coordinates moving with the wave packet at its group velocity ∇ k ω | k . For the dispersion relation of a Boussinesq IGW accord-ing to Eq. (A.9), the group velocity is given by Eq. (A.14) and ∂ k y , k y ω | k = N k , x k ,y − k , x (cid:16) k , x + k ,y (cid:17) / . (A.29) Appendix B: Supplementary plots − − − v r i n c m s − ‘ = 0 f = 921 . µ Hz f = 788 . µ Hz . . . . . r/R ? k s i m / k t h e o Fig. B.1.
Check for the dispersion relation of sound waves. The upperpanel shows the amplitudes of the radial velocity at two di ff erent fre-quencies for (cid:96) =
0. Circles mark radii where our routine detects peaksin the amplitude. The radial wavelengths are then estimated from thedistance of neighboring peaks. For the blue line, the frequency matchesa standing mode and it shows well defined, distinct peaks. The red linecorresponds to a frequency which is between standing modes and showssmall-scale incoherent oscillations with many peaks. In the lower panel the resulting wave numbers from the simulation are compared to the ex-pectation for sound waves k theo = π f / c sound . We find the blue line ingood agreement with theory for radii above the convection zone whilethe red line considerably di ffff