{'fraction_non_alphanumeric': 0.03656649564296767, 'fraction_numerical': 0.015603876537177295, 'mean_word_length': 4.677980209007676, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Total-Variation-Diminishing Implicit-Explicit Runge-Kutta Methods for the Simulation of Double-Diffusive Convection in Astrophysics

We put forward the use of total-variation-diminishing (or more generally, strong stability preserving) implicit-explicit Runge-Kutta methods for the time integration of the equations of motion associated with the semiconvection problem in the simulation of stellar convection. The fully compressible Navier-Stokes equation, augmented by continuity and total energy equations, and an equation of state describing the relation between the thermodynamic quantities, is semi-discretized in space by essentially non-oscillatory schemes and dissipative finite difference methods. It is subsequently integrated in time by Runge-Kutta methods which are constructed such as to preserve the total variation diminishing (or strong stability) property satisfied by the spatial discretization coupled with the forward Euler method. We analyse the stability, accuracy and dissipativity of the time integrators and demonstrate that the most successful methods yield a substantial gain in computational efficiency as compared to classical explicit Runge-Kutta methods.


Introduction
Numerical hydrodynamical simulations are a common tool in astrophysical research. Just as some of their counterparts in the atmospheric sciences and in oceanography, astrophysical fluid flows are characterized by a vast range of timescales which are present in the solutions of the dynamical equations governing the temporal evolution of such flows [23]. Large relative changes of the solutions typically occur on the hydrodynamical timescale τ fluid = ∆x/|u|. Here, u is the local flow velocity and ∆x is the local spatial resolution of the simulation, which coincides with the grid size obtained from spatial discretization of the governing partial differential equations. However, some of the physically important processes can also operate on much shorter timescales than τ fluid . Examples include radiative transfer, sound waves, magnetohydrodynamic processes, and chemical or nuclear reactions (see [23] for example, and references therein).
In stellar astrophysics the two most important among those timescales are that of radiative energy exchange at the scale of a grid cell, τ rad , and the time τ sound a sound wave needs to cross such a cell.
As long as sound waves are energetically or dynamically unimportant, a numerical simulation can be advanced with much larger time steps by using semi-implicit time integration methods, for example, by a fractional step approach (see [11] for a general introduction).
Similarly, if radiative transfer has the numerical characteristics of a stiff problem, as is the case, for instance, for numerical simulations of the surface layers of A-type stars [29], implicit time integration appears desirable as well. Another important example where radiative diffusion can limit the time-step is the numerical simulation of double-diffusive processes in stellar interiors.
Semiconvection is the most important special case of double-diffusive convection in astrophysics. Models of stellar structure and evolution predict settings where the heavier product of nuclear fusion provides stability to a zone which otherwise would be unstable to convective overturning, because temperature sufficiently rapidly decreases against the direction of gravity. Such a zone would become convective if its composition were mixed. The question whether such a zone should be treated as if it were mixed or not is referred to as the semiconvection problem (see [26], [38], and Chap. 13.3 and 13-A in [47], for example). A thorough physical analysis of the semiconvection problem based on numerical simulations in two spatial dimensions for a parameter set relevant to stellar astrophysics is given in [49]. Further discussions and reviews on this topic can be found, for instance, in [4,5,24,43,46]. For this problem, long total integration times are required even when the microphysics is idealised, whereas the time integration step is governed by τ rad and τ sound .
For reasons outlined above, in this paper we discuss the advantages of implicitexplicit (IMEX) Runge-Kutta methods for simulations of stellar convection and diffusion in the parameter regime commonly associated with semiconvection as discussed in [49]. These methods treat only part of the right-hand sides implicitly, where the resulting (generally nonlinear) equations can be solved by means of a generalized Poisson problem. It turns out that the total-variationdiminishing (TVD) property is essential for a numerical time integrator to be successful in simulations of the problems in our focus: to suppress spurious oscillations in the spatial discretization (which in this paper we realize for the hyperbolic terms by essentially non-oscillatory schemes and by dissipative centered finite difference schemes for the parabolic terms), this property has been demonstrated to be necessary for a stable integration in [13].
The TVD property is more generally referred to as strong stability preserving (SSP) or monotonicity when norms other than the total variation norm or even sublinear functionals are considered. When the space discretization has the property that the functional of the discrete spatial profile is decreased in the course of numerical time propagation by the forward Euler method for a time-step ∆t FE , then an SSP method preserves this property under a step-size restriction of the form ∆t ≤ C∆t FE with C > 0. Since the term totalvariation-diminishing is more commonly used in the context of astrophysical simulations, where the total-variation-seminorm is the functional of interest, we mostly use these terms here synonymously. We expect the SSP (or TVD) IMEX methods to be useful also for other astrophysical problems where a high radiative (conductive) diffusivity of internal energy (temperature) restricts the time-step of hydrodynamical simulations, such as simulations of stellar surface convection with steep temperature gradients or at high resolution.
The outline of the paper is as follows. First, we introduce SSP IMEX Runge-Kutta schemes and survey the related literature in Section 1.
Next, in Section 2 we specify the general set of equations to be solved in numerical simulations of semiconvection and related flows and describe the general solution techniques implemented for this class of problems in the ANTARES code [35] which we use for the numerical examples discussed further below. Subsequently, we discuss how this framework has to be modified when solving the dynamical equations with the IMEX approach.
In Section 3 we analyse several SSP IMEX methods from the literature with respect to their radius of absolute monotonicity, stability and dissipativity. We show that the methods yield a significant advantage over classical explicit Runge-Kutta schemes with respect to both efficiency and accuracy and we also suggest a modification for one of the methods, which turns out to improve its efficiency.
We then present numerical simulations of semiconvection in Section 4 to demonstrate the efficiency of the SSP IMEX methods as compared to the classical, explicit SSP Runge-Kutta time integrators and some non-SSP IMEX methods from the literature by giving numerical examples for a single layer in a physical scenario similar to that one studied in [49]. We conclude this paper by a summary of the main properties of the IMEX methods, suggesting reasons for preferring particular methods and providing an outlook on interesting applications, which appear especially suited for this numerical approach.

Implicit-Explicit Runge-Kutta Methods for Semiconvection
To introduce our numerical methods in an abstract setting, we consider the ODE initial value probleṁ y(t) = F (y(t)) + G(y(t)), y(0) = y 0 , where we assume that the vector fields F and G have different stiffness properties. For this type of problems, partitioned Runge-Kutta schemes [15], also called additive Runge-Kutta schemes, are popular. Methods of this kind use different Runge-Kutta formulae for the treatment of the two vector fields. We will see that the spatial semi-discretization of (15) below and the associated boundary conditions give rise to this kind of system.
An s-stage partitioned Runge-Kutta method characterized by coefficient matrices A = (a i,j ) andÃ = (ã i,j ) defines one step y old → y new by If a i,j = 0 for j ≥ i, the method is referred to as an implicit-explicit (IMEX) method.
These methods have first been investigated with respect to the SSP or TVD property in the context of hyperbolic systems with relaxation, where G = 1 εĜ , ε 1, in [36]. Ibidem, the common specification for strong stability preserving IMEX methods is introduced. An IMEX method is referred to as 'SSPk(s, σ, p)' when it has the following properties: k is the order of the method in the stiff limit ε → 0, which is characterized by the coefficients for the explicit part. The latter must necessarily be SSP and is referred to as the asymptotically SSP scheme. s is the number of stages in the implicit scheme and σ the number of stages in the explicit scheme. p is the global order of the resulting combined method. It is essential to observe that if the implicit scheme characterized byÃ = (ã i,j ) is a diagonally implicit Runge-Kutta (DIRK) method, then the explicit part is evaluated only once in each stage, providing the desired computational advantage [36].
The analysis in [36] is valid only for ε 1 [21]. However, several useful examples of strong stability preserving IMEX Runge-Kutta methods are given, see Section 3.
[21] develops a comprehensive theory of strong stability preserving additive Runge-Kutta schemes which extends the concepts for standard Runge-Kutta methods in a natural way: Let τ,τ be the step-size restrictions for monotonicity of the explicit Euler method for the vector fields F and G, respectively. We define the region of absolute monotonicity where the absolute monotonicity at a point (r 0 ,r 0 ) is characterized by algebraic relations for the matrices A,Ã. The boundary in the first quadrant, ∂R(A,Ã)∩{(r,r) : r,r ≥ 0}, is denoted as the curve of absolute monotonicity. The significance of the region R(A,Ã) is expressed in the following theorem [21]: Theorem 1.1 Let (A,Ã) be absolutely monotonous at (−r, −r) with step-size coefficients τ,τ . Then for h ≤ min {rτ,rτ }, diminishing of the norm holds, y i ≤ y old , i = 1, . . . , s, y new ≤ y old .
[22] gives order barriers for strong stability preserving additive Runge-Kutta methods similarly to [28]. The order of an additive Runge-Kutta method (A,Ã) is bounded by the orders of A andÃ, respectively. This implies for IMEX methods the order barrier p ≤ 4 [28]. Moreover, [22] gives a simple algebraic criterion for a nontrivial region of absolute monotonicity in terms of incidence matrices of A,Ã. Some examples of strong stability preserving IMEX Runge-Kutta methods analysed in Section 3 can be seen in [22,28].
Some other relevant issues can be studied for IMEX methods. In this paper we focus on stability regions, error constants and dissipativity analysis. We close this section with a brief description of these concepts.
The stability region of IMEX Runge-Kutta methods is defined in [1,2] via the test equation of the form (1), where For this problem, and the stability region is the part of the complex plane where |R(z)| < 1.
We will perform the corresponding analysis of the stability function for the methods considered in Section 3.
To determine error constants of the methods we have computed the empirical convergence orders by solving the non-linear test problem y (t) = (1 + sin(y(t))) + (y 2 (t) − sin(y(t))), with the known exact solution y(t) = tan(t). In this paper, the error constants are determined from the errors at t = 1.3. Their size is vital for the comparison of the accuracy of the methods of the same order and therefore the assessment of the work/precision relation. Of course, this single example only gives a rough indication of the size of the error constant and no rigorous estimate, but it seems sufficient for our purpose of comparing the methods in our focus with respect to accuracy, which we do in Sections 4 and 5 further below. The example was chosen such as to represent a nonlinear initial value problem with known solution whose profile is arbitrarily unsmooth as t → π 2 .
Finally, we study the dissipativity of the time integrators in conjunction with suitable space discretizations. [45] gives a justification for considering only the diffusion term in this context since the advection term becomes negligible asymptotically. We will thus investigate the dissipativity of the implicit scheme specified byÃ. To this end, we apply the spatial discretizations L ∆x in our focus to the heat equation u t = bu xx , and associate for the spatial discretization u j±k ↔ e ±iθk . Thus, we compute the amplification factor g(θ, µ) = R((L ∆x u) j ), with µ = b ∆t (∆x) 2 . This represents the factor by which oscillations of frequency θ are amplified in each time step. We pay particular attention to the case θ = π corresponding to the mesh width. If g = 0 or |g| = 1 also for a smaller 0 < |θ 0 | < π, then this value would represent the limit for a robust integration. However, such a pathological behaviour is only conceivable for methods of higher order.
The spatial discretizations which were found to show a dissipative behaviour in [27] are the second-order three-point scheme (the upper index refers to the time step) and the fourth-order stencil These are the methods actually implemented in ANTARES, where (9) is the default (see also [27]).

Solving the Hydrodynamical Equations with ANTARES
In our model, the fundamental equation of motion is the fully compressible Navier-Stokes equation which describes momentum conservation: The state variables in the model equations generally depend on the spatial variables (x, y, z) and time t. In the simulations presented in Section 4 we solve problems in two spatial variables, whence the variable z will be dropped in the rest of the paper. The (explicit) dependencies are stated in Table 1. For simplicity, we omit the dependencies in the problem specification (10)- (14). The model is completed by the continuity equation which ensures conservation of mass, and the total energy equation e + ∇ · (u(e + p)) = ρ(g · u) + ∇ · (u · σ) + Q rad , which describes conservation of the latter. In the case of a two-component fluid, the system is augmented by the concentration equation of the second species, (cρ) + ∇ · (cρu) = ∇ · (ρκ c ∇c).
The variables and parameters which appear in the model formulation are collected in Table 1.
In general, the radiative source term Q rad is determined as the stationary limit of the radiative transfer equation which is solved for all ray directions r and for all frequencies ν, resulting in the specific intensity I ν (r), for details see [47]. S ν here denotes the source function.
The equations of hydrodynamics (10), (11) and (12) are closed by the equation of state which describes the relation between the thermodynamic quantities.
For the particular choice, see [35].
For the initial condition, a slightly perturbed static model stellar atmosphere or stellar envelope is used which is equipped with a small seed velocity field or density perturbation to start dynamics away from equilibrium.
In the framework discussed below, boundary conditions are based on the assumption that all quantities are periodic in both horizontal directions. Moreover, for the hydrodynamical equations, "closed" (Dirichlet) boundary conditions at the upper and lower boundary of the computational domain are used. A recent development is to replace these by "open" (Robin) boundary conditions. These allow inflow and outflow of fluid along the vertical direction which is defined by the direction of g. For the radiative transfer equation (14), incoming radiation at the boundary of the computational domain must be specified.
Since double-diffusive convection in stars takes place in regions which are optically thick, the quantity Q rad can accurately be obtained by means of the diffusion approximation for radiative transfer, Q rad = ∇ · F rad = ∇ · (K∇T ).
In this case, further knowledge about the intensity I ν is not necessary.
The ANTARES code [35] solves this system of equations numerically in either one, two, or three spatial dimensions on a rectangular grid (spherical coordinates with a logarithmically rectangular grid are also possible, i.e., the grid may be locally rectangular with logarithmic grading in the radial component).
For the spatial discretization, ANTARES allows the definition of several grids which can be nested inside each other to improve resolution in regions of interest. At the moment, ANTARES provides up to three levels of nested grids. For the hyperbolic terms, discretizations of ENO (essentially non-oscillatory [40]) type are implemented. These comprise classical ENO methods, WENO (weighted essentially non-oscillatory) methods [40] (optionally in conjunction with Marquina flux splitting [8]) and CNO (convex non-oscillatory) schemes [31]. Each of the methods uses adaptive stencils which are chosen such as to avoid spurious oscillations in the computed solution. The spatial derivatives are calculated for each direction separately.
The parabolic terms are discretized by dissipative finite difference schemes [27] of fourth order. The radiative heating rate is determined by the short characteristics method, or by means of the diffusion approximation for radiative transfer, where appropriate. For the time integration, total variation diminishing Runge-Kutta methods [39,41] are employed.
ANTARES implements two different parallelization concepts. For architectures with distributed memory, domain decomposition is used and realized by an MPI implementation. In this approach each grid is split along the horizontal direction(s) and optionally, also along vertical ones, into subdomains. The memory required to store the computational variables for each subdomain is provided by the resources available to the CPU core performing the computations necessary for that subdomain. In this way, each CPU core is mapped to a specific geometrical volume. However, since some supercomputers offer only a limited amount of memory per CPU core and because the domain decomposition approach is not very efficient on small grids, ANTARES offers a second type of parallelization which can be used along with or independently of the former. It is based on a shared memory concept for each subdomain and is implemented through OpenMP directives. Therefore, the most time consuming operations which can also be performed independently of each other are identified and computed in parallel. This approach scales only to a moderate number of CPU cores, but allows improvement of the scaling and the computational speed of the domain decomposition based parallelization for a larger number of problems and for a greater variety of computer architectures.
In the following, the dynamical evolution of the fluid is described by the multispecies Navier-Stokes equations presented above. Additionally, dimensionless quantities such as the Prandtl number Pr = c p η/K, the Lewis number Le = c p ρκ c /K, the Rayleigh number Ra and the stability parameter R ρ are defined to determine the diffusivities κ T , κ c and the viscosity η. The former quantities arise in the definition of the starting model but do not appear in the evolution equations (15) below. Since we solve the dynamical equations for a compressible flow, we specify the vertical extent of the simulation domain in multiples of the pressure scale height H p = P/(ρg). For the simulations presented in Section 4 the domain always covers 1H p . In the physical model, intermolecular forces are neglected, so the fluid is assumed to be an ideal gas. The radiative source term Q rad is modelled using the diffusion approximation with a heat conductivity K which is constant in time and otherwise only a function of the vertical coordinate [33,34]. We use a variant of this setting here, where not only Pr, Le, and c p , but also K, κ c , and η/ρ take constant values [49]. This setup simplifies studies of the basic physics while it is still useful for extrapolations to astrophysically relevant cases (cf. [49]).
For the model problem, the multispecies Navier-Stokes equations can be recast as . (15) In the context of the problem (15), the i th implicit stage of an IMEX Runge-Kutta method is typically of the form where y * is known from previous stages. This is a consequence ofã i,j = 0 for j > i.

This translates to
Rearranging (18) leads to Obviously, this is a general elliptic equation for c i of the form Due to model assumptions, (20) can also be transformed to resemble a general elliptic equation. We start by recalling e = e int + e kin (23) Bearing in mind equation (19), equation (20) reads The equation of state for an ideal gas 4 relates the temperature T and the internal energy e int via if we assume the ratio of the specific heats at constant pressure and volume to equal 5/3. Here, m denotes the mean molecular weight of the compound.
So we have at stage i and therefore, we arrive at Since m i is evaluated using the mass fraction c i , it is necessary to solve equation (21) first.
Thus, the solution of an implicit stage translates to the solution of the generalized Poisson equations for the mass fraction c and the temperature T . In the ANTARES framework, finite elements are used for the discretization of (22). The resulting linear system is solved by the conjugate gradient method. For parallel computations, the Schur complement algorithm is applied. A detailed description is given in [14].
The above procedure applies without modification to the case of the fully compressible Navier-Stokes equation. However, for low Mach number flows a splitting approach is preferable where the terms containing pressure are treated separately in a post-processing step, i. e. after evaluation of all other terms for the computation of the velocity fields. The latter are obtained for the current step from an additional generalized Poisson equation for the pressure. This procedure is described in detail in [11]. Consequently, the explicit stage is evaluated here using a fractional step method [30] as implemented in [17].

Strong Stability Preserving IMEX Schemes from the Literature
In this section we study different SSP IMEX methods from the literature focusing on the topics described in Section 1. The results are summarized in Figure 16 and Table 10.
All the strong stability preserving IMEX schemes listed in the following subsections have DIRK (diagonally implicit Runge-Kutta) methods as the implicit scheme. This structure ensures that the stages can be solved successively, and the explicit part only has to be evaluated once in each stage.
The stability region is entirely located in the left half plane, tangent to the imaginary axis and unbounded as (z) → −∞. Hence, the schemes are A(α)stable with α = π 2 , but not A-stable [16]. Moreover, lim (z)→−∞ R(z) = 0. A plot of the stability region is given in Figure 1 (left). The stability function of the implicit schemeÃ is A plot of the related stability region is shown in Figure 1 (right). The schemẽ A appears to be A-stable and satisfies lim (z)→−∞ RÃ(z) = 0, implying Lstability.
The first positive zero of g(π, µ) is ≈ 0.6035 for the three-point scheme (8), where the function changes its sign, and |g(π, µ)| never exceeds 1. The first positive zero of g(π, µ) for the fourth order space discretization (9) is ≈ 0.4526, where the function changes its sign. The modulus never exceeds 1.

Modification of γ
We may conceive of optimizing the method (29) by adapting the value of the parameter γ in the definition ofÃ according to the resulting stability, accuracy, and dissipativity properties. The region of absolute monotonicity depends on γ as follows [19]: A plot of the function R(Ã) in dependence of γ is given in Figure 2.  It was demonstrated with a Matlab implementation that the convergence order two is retained also for the modified values of γ. The error constant depends on γ, however. In Figure 4 (right) we plot the error constant as a function of γ, where the error is determined at t = 1.3 for the test problem (7). We note that for small γ, the error constant decreases as γ grows, while for γ > 0.1833 the constant grows monotonically. This behaviour does not appear to be related to the results we obtain for the dissipativity analysis, where for γ = 0.25 the behaviour changes.  To assess the dissipativity properties of the modified scheme we vary γ and compute the first positive zero of g(π, µ) and the point µ where the modulus of this function exceeds 1.
To illustrate this analysis, in Figure 5 we vary γ and compute for 30 values of γ spaced equidistantly in the interval [0.15, 1.1 − 1/ √ 2] the first positive zero of g(π, µ) and the point µ where the modulus of this function exceeds 1. This analysis is given for the three point space discretization (8) in Figure 5. We conclude that it may be of interest to choose γ < 0.25 to ensure 0 < g(π, µ) < 1. This can be achieved by a value of γ just slightly smaller than 0.25.
Taking into account the results displayed in Fig. 4 we suggest to optimize γ by taking it as large as necessary to avoid any linear stability restrictions due to the term G(y(t)) in (1) and as small as possible to minimize the stability constant C. Note that some special values for γ are 0, (1 − 1/ √ 3)/2 ≈ 0.2113, and 1/4, in addition to the originally proposed value of 1 − 1/ √ 2 ≈ 0.2929. They are readily identified to yield the classical explicit SSPRK(2,2) method 5 of order two by Shu & Osher [41], the optimal implicit third order SSP method with two stages [10], and the optimal implicit second order SSP method with two stages [10]. In our numerical tests reported in Section 4 below, we found that the choice γ = 0.24 yielded the most efficient time integrator.

An SSP2(3,3,2) Method
This is a modification of a scheme from [36], where the latter turned out to have a trivial region of absolute monotonicity. It holds that R(A) = 2 and We note that the latter is a correction with respect to [21], since we have found r to be necessarily bound by 1 in R(A,Ã). A plot of R(A,Ã) is given in Figure 6. A plot of the stability region is given in Figure 7 (left). We observe that the stability region is tangent to the imaginary axis and appears to be unbounded implicit scheme is given by A plot is shown in Figure 7 (right). Note that the stability region is not connected. The method is A-stable, however. Moreover, the schemeÃ satisfies lim (z)→−∞ RÃ(z) = 0.
The dissipativity analysis for the implicit scheme defined byÃ yields the amplification factors for the standard three-point space discretization (8) and the fourth order stencil (9). These amplification factors are evaluated at the points θ ∈ {0, π 4 , π 2 , π} in Tables 4 and 5, respectively.

An SSP3(3,3,3) Method
The stability region occupies a bounded domain in the negative half-plane, and slightly overlaps the imaginary axis, see Figure 9 (left). Note that A plot is shown in Figure 9 (right). The point where the stability region intersects the negative real half-line is located at x ≈ −3.248 for the IMEX scheme. The same value is computed for the stability region of the implicit scheme. However, the stability region ofÃ extends further along the imaginary axis. Finally, lim (z)→−∞ |RÃ(z)| = ∞.
The dissipativity analysis for the implicit scheme defined byÃ yields the amplification factors for the standard three-point space discretization (8) and for the fourth order stencil (9). These amplification factors are evaluated at the points θ ∈ {0, π 4 , π 2 , π} in Tables 6 and 7, respectively.

Numerical Experiments
In this section, we present the results of the simulations performed with the time integrators discussed in this paper and compare their performance to the classical explicit 2-stages second order and 3-stages third order methods of Shu & Osher [41], and the 3-stages second order method of [28]. They are the optimum SSP RK methods for their given number of stages and order and we refer to them by their common technical specifications 'SSPRK(2,2)' and 'SSPRK(3,3)', as well as 'SSPRK(3,2)', respectively. The coefficients of the SSPRK(2,2) and SSPRK(3,3) methods were derived for a different purpose in [18] and in [9], where the third order method was proposed as an embedding formula for the second order method (see also [3]). Shu & Osher [41] derived a first framework for deducing higher order Runge-Kutta methods with the total variation diminishing property and first identified the optimal explicit second and third order methods with two and three stages, respectively. For these SSPRK(2,2) and SSPRK (3,3) schemes an analysis of their stability, dissipativity and accuracy properties can be found in [27]. For the SSPRK(3,2) method, a similar study is given in the Appendix section further below. Note that SSPRK(2,2), SSPRK (3,2), and SSPRK (3,3) are the explicit schemes in the IMEX methods (29), (33), and (35), respectively.
Furthermore, we also consider some non-SSP methods from the literature, both explicit and IMEX, for the sake of comparison, and also an asymptotically stable SSP IMEX scheme from [36].

Implementation Issues
To define the test problem, according to [33] we specify a hydrostatic configuration which is unstable against convection. The simulation of a single semiconvective layer requires the mean molecular weight to be linearly (and stably) stratified. As time evolves, we expect convection to set in and mix the zone completely, although its development is inhibited by the stable mean molecular weight gradient. A critical quantity in this process is hence the buoyancy timescale and we return to this topic further below.
The simulations shown here have been performed on the Vienna Scientific Cluster, using 64 CPU cores in parallel. The spatial resolution is 400×400 grid points. Simulation time is measured in units of sound crossing times (scrt). where with (not necessarily equal) Courant numbers (CFL numbers) C c , C T , C visc , C fluid . This assumes that the time-step limitation due to sound waves has been removed by a fractional step approach as mentioned at the end of Section 2. Otherwise, τ fluid additionally depends on the sound speed. We note that the source term in F (y(t)) in (15), which represents buoyancy forces acting on the flow, can be neglected in the limit where max {∆x, ∆y} → 0, since its contributions are of lower order (see [45] for a discussion of the treatment of lower order terms in stability analyses).
Due to (15), IMEX methods treat the terms ∇ · (ρκ c ∇c) and ∇ · (K∇T ) implicitly, so the restrictions τ c and τ T do not have to hold. Since at least the first part of simulations of semiconvection is usually dominated by diffusion processes and the Prandtl and Lewis numbers satisfy Pr < 1, Le < 1 in a stellar context, the simulations are initially restricted by τ T . Hence, IMEX methods can provide the desired computational advantage.
To enhance the stability and the efficiency of our methods we have implemented a heuristic to adaptively select the time steps. The most effective criterion for regulating the time steps turned out to be to monitor two-point instabilities appearing in the conservative variables (ρ, ρc, ρ u, Et). Due to the gravitational force operating vertically, such instabilities are prone to appear in the horizontal direction.
To detect the occurrence of such oscillations, for each variable we use the difference between two grid cells to determine the sign of the corresponding gradient. In case of the density ρ this reads where 1 < i < n x and 1 < j < n y , assuming the grid consists of n x ×n y points. If the sign pattern of (d 1 , d 2 , d 3 , d 4 ) corresponds to (+, −, +, ±), (−, +, −, ±), (±, +, −, +, ) or (±, −, +, −), we have located a two-point instability. Since the time-step is initially chosen to be that one required for a fully explicit time integration method as in (37), such patterns are smoothed out rapidly, if present in the initial condition. Consequently, their later occurrence is a good indicator for an instability developing because of too large a time-step taken during the time integration.
The time-step control permits the occurrence of n y · 0.1 two-point instabilities for fixed i. If this limit is exceeded, the time-step is repeated using a step-size decreased by a factor 2 3 .
To permit the system to readjust after reducing the time-step no modifications of ∆t are made for the next 15 time-steps regardless of the number of twopoint instabilities. If the number of oscillations still exceeds the given limit after those 15 time steps, the time-step is again reduced. If no or very few two-point oscillations are encountered for over 50 successive time-steps, the time-step is augmented by a factor of 5 4 .
Alternatively to this heuristic control, we could also monitor the rate of change in the solution to adjust the time-steps. However, the proper rate required to prevent the development of two-point instabilities for the present application turned out to be too pessimistic to achieve CFL numbers as high as discussed below. Furthermore, a simple control of the residual did not yield a satisfactory  Table 8 Comparisons of time steps and CFL-numbers over the first 80 scrt for the case where Pr = 0.1 and over the entire 200 scrt for the case where Pr = 0.5 (see also text).
behaviour. This leaves the heuristic time-step control as the most effective method for the IMEX based time integration of our simulations. Tables 8 and 9 sum up the performance achieved with the presented Runge-Kutta schemes. The evolution of the time steps is illustrated graphically in Figures 10 and 11 for the two simulation scenarios we have focused on. Since the capability of the method is best judged in that part of the simulation where the fluid velocity is too small to severely limit the time-step ∆t, the CFL-numbers and time-steps listed in Table 8 have been measured in this regime. For the first test scenario this corresponds only to the first 80 scrt (note the slope beginning after about 100 scrt in Figure 10).

Numerical Results with SSP Schemes
In Table 8 we compare the largest possible time steps ∆t max , the average time steps ∆t mean , the maximal and average CFL numbers resulting from the adaptive step selection, and the initial CFL numbers. The tests were performed for Prandtl numbers Pr = 0.1 and Pr = 0.5 distinguishing Simulations 1 and 2, respectively, a Lewis number Le = 0.1, R ρ = 1.1, and a modified Rayleigh number Ra * = 160000 related to the Rayleigh number through Ra * = Ra · Pr.
Comparing the performance of the second order schemes it is obvious that the modification of γ in IMEX SSP2(2,2,2) has a stunning effect on the stability of the scheme. The positive definiteness of dissipation in this method proves most effective in suppressing oscillations, permitting an average time-step and CFL-number up to two thirds higher than the original IMEX SSP2(2,2,2) method. Table 8 shows that IMEX SSP2(3,3,2) permits a time-step and CFL number more than twice as high as IMEX SSP2(2,2,2) even with modified γ. However, a comparison of the computation time given in Table 9 shows that this does not improve by the same factor as the CFL-number, since as the time-step ∆t grows, the iterative solver for the generalized elliptic problems requires more iterations to converge, resulting in an increase in computation time. A comparison of the computation times shows that, although IMEX SSP2(3,3,2) permits impressively large time steps, the need to solve three additional generalized Poisson problems related to the third stage takes its toll, whence the method's performance is inferior to IMEX SSP2(2,2,2) with modified γ and in case of Simulation 2, is not even competitive to SSPRK (3,2), though it still performs better than the SSPRK(2,2) scheme. Note that for Pr = 0.1, the IMEX methods outperform even the best explicit method, while for the more moderate Pr = 0.5, the best explicit integrator SSPRK(3,2) is slightly more efficient whereas the classical methods noticeably lag behind.
Interestingly, the initial (preset) CFL number has a negligible influence on the actual CFL number reached in the diffusive part of the simulation. However, as soon as the fluid velocity seriously restricts ∆t, a higher initial CFL number leads to a significantly larger average time-step and reduces the required  Table 9 Comparisons: computation times and overall number of time steps over 200 scrt. computation time, since its value is used to define the time-step restriction for the terms integrated with the explicit part of the IMEX scheme.
Once the time-step is limited by τ fluid , the remaining gain in ∆t by the IMEX schemes in ANTARES is essentially due to the optimization of the time-steps by the algorithm explained further above. Since ∆t is rather small in that case  the convergence of the generalized Poisson solver is fast enough to allow the IMEX schemes to lag only slightly behind their explicit counterparts.
We point out that optimization of the solver for the generalized Poisson equation (21) has the potential for a further significant decrease of computation time required, both in the diffusive regime, but also if the time-step is limited by τ fluid . This is important since changing the time integration method during a simulation run is not advisable if both the diffusive and the convective phase should be interpreted in a consistent manner.
A comparison of SSPRK (3,3) and IMEX SSP3 (3,3,3) also shows that the larger time-steps of the semi-implicit method lead to a significant gain in computational efficiency in case of a third order method.
We have investigated the accuracy of the time integration with IMEX schemes by a comparison to a reference solution obtained with the SSPRK(2,2) method. The reference solution was computed on the same spatial grid but with a time-step eight times smaller than that one mentioned in Tables 8  and 9 for this method, i.e. for a CFL-number of 0.025. Figure 12 shows the root mean square difference of the mass ratio He / (H + He), obtained by summation over all grid points and normalization relative to their number, between the numerical solutions of the second order SSP IMEX methods and the reference solution for each case. For the IMEX SSP2(2,2,2) method both the results for the standard choice of γ = 1 − 1/ √ 2 and the best performing value of 0.24 are displayed. We also show the normalized root mean square differences between the reference solution and the SSP second order explicit methods computed with their standard CFL number given in Table 8.
For both Simulation 1 and 2 one can easily spot the initial increase of the error due to the growing time step for the IMEX methods induced by the automatic time step control. A plateau is reached once the time-step stabilizes around a typical mean value (cf. also Figures 10 and 11). Note that the IMEX SSP2(2,2,2) method with γ = 0.24 has a smaller error than the original IMEX SSP2(2,2,2) method, as expected from the error constant shown in Figure 4. The largest differences occur for the IMEX SSP2(3,3,2) method which also has the largest mean time-step. The error constants of the different methods (see also the summarizing Table 10 below) provide a rough measure for comparing simulation runs with similar time-steps.
However, a comparison for a given point in time has only limited meaning. One of the reasons is that the solution changes its nature as a function of time. Initial vertical oscillations are damped out (at least first few scrt), then the velocity field slowly starts building up (visible after 15 scrt), followed by the formation of large scale gravity waves (oscillatory behaviour of the error in the range between 25 and 100 scrt) until the waves start to break   To further illustrate that the large time steps of the IMEX methods during the diffusive phase do not degrade the accuracy of the time development of the solution, we compare the simulation results obtained by the different time integrators in Figures 13 and 14. The pictures show the mass ratio of He vs. He+H at each spatial point at a given instant in time. Figure 13 demonstrates that after the end of the diffusive phase, just at the onset of turbulence, which occurs at ∼ 125 scrt for this problem, the results are quite comparable. The most visible differences can be found for the case of the simulation with the largest time-steps (picture (d)). It also shows that the root mean square errors displayed in Figure 12 are negligible on a qualitative (and rough quantitative) level as long as they are smaller than about 10 −3 . As expected, however, some time after the onset of turbulence the solutions necessarily have drifted apart. This is demonstrated in Figure 14, where the solutions already look different from each other.
Recalling Figure 12 it is not surprising to find the largest differences in Figure 13 for the schemes having the largest time-steps during the diffusive phase. Still, the large scale structures of the solution begin to diverge only once the turbulent phase has been reached which in turn for each of the different time integration methods occurs after about the same integration time t. Large time-steps during the diffusive phase hence yield acceptable accuracy. Indeed, the spatial resolution is more important than the temporal one. This can be demonstrated by using a high resolution grid of 800×800 points. We have performed such a reference run for the case of Simulation 1 with the SSPRK(2,2) method for time integration. Note that the doubling of resolution leads to a four times smaller time-step during the diffusive phase. Looking at the first row (pictures (a) and (b) of Figure 15) the differences at 100 scrt, i.e. just after the onset of wave-breaking and the beginning of the turbulent phase, are still small: in the simulation with higher spatial resolution the breaking tips are somewhat more pronounced and the contrasts sharper. This has changed at 125 scrt shown in the second row of Figure 15. The common initial condition may still be inferred, but the simulations have notably evolved away from each other. Clearly, the spatial resolution is much more important than the influence of the time steps and the time integration method chosen, since picture (c) of Figure 15 is essentially indistinguishable from its counterpart with standard time resolution, picture (a) of Figure 13. Furthermore, at 100 scrt the simulation using the IMEX2(3,3,2) method for time integration, which has the largest time-step during the diffusive phase, is nearly indistinguish-  able from the SSPRK(2,2) run with high time resolution shown in picture (a) of Figure 15 (the IMEX results are not shown here for the case of 100 scrt, since the differences to the latter picture are very difficult to spot).
We conclude that the spatial resolution is indeed more important than high temporal resolution and large time steps during the diffusive phase are clearly tolerable for simulations of astrophysical convective flows, if they do not affect stability. Resolutions of 800 grid points per spatial direction in 3D are usually not affordable anyway and sometimes even large parametric studies in 2D may still be too expensive (for instance, for the case of semi-convection and purely explicit time integration methods). We note that the necessity of a resolution of 400 points, which has been used for most of the simulation runs shown here, was calculated following [48] and [49] where in turn the physical arguments of [43] had been used to estimate the thickness of solute and thermal boundary layers in semi-convection and the applicability of this approach to the parameter range we are interested in had been confirmed. Thus, for the parameters of the more demanding one of our models, Simulation 1, we concluded the smallest structures of interest, the solutal boundary layers, to span 6 grid points, if the whole box is discretized by 400 points in each direction. Examples for them are the top and bottom boundary layers which can easily be seen in Figures 14 and 15. Indeed, in the simulation with a high spatial resolution of 800 points in each direction these layers are hardly any thinner than in the case of 400 points (cf. the bottom row of Figure 15 which shows the simulations developed well into their final, turbulent state).
The duration of the diffusive phase is determined by the stability of the stratification, parametrized through R ρ , and related to the buoyancy term, the second term on the right-hand side of (15) and thus also the second term of F (y(t)) in the same equation. Timescales related to this term can be computed for a variety of physical problems. They include the reciprocal of the growth rate of small density perturbations when a fluid of higher density ρ 2 is layered above fluid of lower density ρ 1 (this situation is denoted as Rayleigh-Taylor instability, see [6]). In this case the growth rate follows the dispersion relation ω 2 = gk(ρ 2 − ρ 1 )/(ρ 2 + ρ 1 ) for a local gravitational acceleration g and a perturbation with wave number k. Its magnitude can be bounded by the simple relation ω 2 = gk. Similar dispersion relations are found for gravity waves and growth rates of convective instability (see the classical paper [7] and also [26] for a summary). As implied by ω 2 = gk and the form of (15), one can estimate a buoyancy time-step restriction by t buoy = min{(∆x) 1/2 }/g 1/2 (see also [32]). For the present simulations we find t buoy ∼ 360 s, i.e. about 0.069 scrt.
Note that this is about four times larger than the largest time-step reported in Table 8. Though irrelevant in the asymptotic limit and not a constraining quantity here, one might still consider this term to be integrated implicitly in other cases. However, if the excitation and breaking of waves is important to describe the onset of the turbulent convective flow, as in the present case, damping of such waves by an implicit time integration may be undesired. Hence, explicit time integration of the buoyancy term could be preferred for physical reasons, even if the time-step were actually constrained by such a splitting for the time integration.

Numerical Results with non-SSP Schemes
From Tables 8 and 9 one can readily see that during the diffusion dominated phase the SSP IMEX schemes achieve time-steps which exceed the region of absolute monotonicity ensured by (32) for IMEX SSP2(2,2,2) and their counterparts given after (33) and (35) for IMEX SSP2(3,3,2) and IMEX SSP3 (3,3,3), respectively. One might hence question whether the property of absolute monotonicity is really necessary for the time integration of the numerical simulations we have considered above. To show that this property is indeed required we have performed several test runs for the case of Simulation 1 with time integration schemes which do not diminish the total-variation norm.
The first candidate we have investigated is the ARK3(2)4L [2]SA scheme proposed in [25]. This is an L-stable, stiffly accurate third order, additive Runge-Kutta method with four stages. With its choice of coefficients it belongs to the group of IMEX methods. However, neither its explicit nor its implicit part are strong-stability-preserving (the Butcher arrays feature negative coefficients, cf. Theorem 4.2 in [28]), hence also the entire method does not fulfill the criteria of Theorem 1.1 on absolute monotonicity. If the SSP-property were of no importance, this method should be quite robust. However, it turns out that this is not the case when we use it to integrate Simulation 1 in time. Taking CFL start to be 0.2 as for the other IMEX methods presented in Table 9, the time integration with ARK3(2)4L[2]SA crashes after just 78 time-steps. Indeed, CFL start has to be lowered to 0.1 to successfully launch the simulation. But over the first 10 scrt CFL max is found to never exceed ∼ 0.2 and the average CFL mean is only ∼ 0.15. In conclusion the size of the time-steps achieved with this method have been found to not exceed those achieved with the explicit, three-stage, third-order SSPRK(3,3) method of [41]. Compared to IMEX SSP2(3,3,2), the maximum and mean CFL numbers are 25 and 20 times smaller, respectively. We have thus given up this simulation run after 10 scrt: evidently, the ARK3(2)4L[2]SA method is not efficient for the kind of problems we are interested in. The implicit, L-stable and stiffly accurate nature of this scheme is not sufficient to provide any advantages on its own during the diffusive phase of the simulation.
To further investigate the importance of the strong-stability-preserving property we selected the classical, explicit, third order, three-stage Runge-Kutta method first proposed in [18] and known as Heun's third order method. This is a non-SSP scheme since not all of its stages are used in the final integration which yields y new , as pointed out in [28], where it was used to illustrate the growth of solutions measured in standard norms for both parabolic and hyperbolic problems in situations where the exact solution is not growing in these norms. By comparison the SSPRK(3,3) scheme was found to not exhibit such growth. When we apply Heun's third order scheme to integrate Simulation 1, we achieve a stable simulation over the entire extent of 200 scrt with the same average time-step and with the same Courant number as for the SSPRK(3,3) scheme. However, during the diffusive phase in Simulation 1 the solution itself is slowly growing in time while a velocity field is being built up by the convective instability and at the grid scale the dissipation is provided by the parabolic terms during the entire simulation.
The semi-convection problem discussed here is a rather benign example for numerical simulations in astrophysical applications. If we apply the same scheme to a simulation of solar surface convection as in [35], i.e. for a case of 219 × 159 grid points and a standard, moderately low resolution of 18.57 × 40 km 2 , and a standard choice for the microphysics with non-grey radiative transfer for the calculation of Q rad , the differences in stability become apparent. For this physical problem the term representing viscous dissipation in the momentum equation does not provide sufficient dissipation on the grid-scale at any resolution achievable in the foreseeable future and the dissipation properties of the temporal and spatial discretization of the advection operator become important (the term implicit large eddy simulation is used in such cases). While the SSPRK(3,3) method and also the SSPRK(2,2) method, used together with the spatial discretization of [40], have no problems in completing a simulation of 20 scrt with a CFL number of 0.25, Heun's third order method leads to a crash after 9.2 scrt. 6 Such failures are usually caused by numerically induced fluctuations in the low density and temperature region of the simulation box which result in negative or at least unphysically small values, whence they fall outside the tabulated region of microphysical properties. We consider this finding as sufficient to exclude non-SSP methods from being recommendable for numerical simulations of stellar convection, since also here we have chosen a rather benign test case within its class. Numerical simulations of stellar surface convection in white dwarfs, A-type stars, or Cepheids reach far more extreme conditions with up to three times higher, super-sonic Mach numbers, density contrasts around shock fronts higher by an order of magnitude and more, and for the case of A-stars and Cepheids, at lower effective resolution because of four to ten times steeper gradients and limited computational resources.
We finish our study of non-SSP Runge-Kutta schemes with a third case: an IMEX Runge-Kutta method where both the explicit and the implicit part are strong-stability-preserving, but the combined scheme does not fulfill the criteria for absolute monotonicity in the sense of Theorem 1.1. The scheme was proposed in [36] (Table IV, page 139) and indeed the IMEX SSP2(3,3,2) scheme (33) is a modification of that scheme proposed in [21] to obtain a nontrivial region of absolute monotonicity. The original scheme of [36] differs from that one only by having the entries {1/4, 0, 0} and {0, 1/4, 0} in the first two rows of the Butcher array of the implicit scheme instead of {1/5, 0, 0} and {1/10, 1/5, 0}, respectively. This scheme performs quite well. Indeed, during the diffusion dominated phase of Simulation 1 its time-steps are even 6% to 8% larger than those achieved with IMEX SSP2(3,3,2), while during the turbulent phase they fall back to at most the size achieved by the scheme (33). Once more, however, we recall that astrophysical simulations often have to deal with limited resolution at least in part of the simulation domain. To reproduce such a case we have run Simulation 1 with both (33) and the original scheme of [36] for the case of only 100 × 100 grid points while leaving everything else unchanged. In this case the boundary layer due to concentration c is represented vertically only by one to two grid points (see the discussion on resolution and reference solutions further above). While during the diffusive phase no major differences become apparent, the behaviour found for the advection-dominated, turbulent phase was discrepant: whereas nothing suspicious occurred for the scheme (33), the time-step of the scheme of [36] dropped to arbitrarily small values after ∼ 113 scrt, which indicates the occurrence of two-point instabilities, and the simulation had to be terminated.
We conclude that only Runge-Kutta methods which are strong-stabilitypreserving have the necessary prerequisites for stable time integrations of astrophysical convection simulations. If, in addition, the time-step is limited by diffusion processes, this limitation can be overcome by IMEX methods provided their explicit and implicit parts are strong-stability-preserving. To ensure stability also in cases of low resolution IMEX SSP methods should also have a nontrivial region of absolute monotonicity as required by Theorem 1.1.

Conclusions and Outlook
In this paper we have given an extensive discussion of the mathematical properties and practical usefulness of total-variation-diminishing implicit-explicit Runge-Kutta methods for the time integration of advection-diffusion equations arising in the simulation of double-diffusive convection in astrophysics.
In this section, we summarize the results obtained in Sections 3 and 4 (stability, dissipativity, accuracy and efficiency), and give a brief outlook on future developments.
The stability regions for the IMEX methods, their implicit sub-parts and the explicit schemes are given for comparison in Figure 16. The left boundaries of the stability regions z left , the points where the amplification factors from the dissipativity analysis become zero and their moduli exceed 1, and the error constants C for all the methods we have investigated are summed up in Table 10.
We found that methods introduced in [20,21,22,36] excel over the classical explicit methods [41]. It was found that among explicit schemes, only the explicit SSPRK(3,2) scheme first proposed in [28] is competitive in situations where explicit time integration can be expected to yield sufficient efficiency and accuracy. Examples for such a scenario include simulations of solar granulation at moderately high spatial resolution, where the time-step limitation associated with diffusion is negligible (see [35] for results on this problem obtained with the ANTARES code and [44] for a review on the underlying physics). From Figure 16 and Table 10, clearly there is no single scheme which features the most advantageous properties in all considered aspects. However, we found in numerical experiments that the most efficient method seems to be (29) with the choice γ = 0.24. This value deviates from the optimal value for strong stability, but leads to a scheme with favourable dissipativity, stability, and accuracy properties. Depending on the domain of stability required for a given problem the value of γ in (29) may be optimized such that it is sufficiently large for stability, but small enough to minimize the error constant, while showing favourable dissipation properties (a strictly positive amplification factor with modulus less than 1 for any wave number k other than zero).
We note that for numerical problems arising from a method of lines approach to the equations of hydrodynamics, as discussed in this paper, lower order methods usually have sufficient efficiency to be competitive, since the spatial discretization limits the overall accuracy. Hence, the best explicit scheme we have tested for this kind of application is SSPRK (3,2), as it permits the largest CFL numbers among methods of this class at an affordable computational cost and with sufficient accuracy. By comparison, the classical methods of second and third order [41] offer the convenience of being usable together as an embedding formula. However, this approach is more than twice as expensive as SSPRK (3,2), as can be seen from a comparison with SSPRK(3,3) using Table 9.  Table 10 Summary of the analysis of SSP integrators. The asterisk in the third column indicates a change of sign at µ for g 4th (±π, µ). Other details are given in the text.
We have also demonstrated that the larger time-steps achieved by SSP IMEX methods reduce the accuracy of the solution during the diffusive phase of the semi-convection simulations by an acceptably small amount. For the applications shown here, and indeed for a majority of astrophysical fluid dynamical simulations, the accuracy is limited by spatial resolution (and thus eventually by existing computational resources) while the time-steps are limited by stability. This makes IMEX methods attractive, since quite often the most severe limitations stem from stiff terms representing diffusion processes (for restrictions due to sound waves other operator splitting based methods are existing). However, as we have shown by a comparison with results from non-SSP methods, it is important that the IMEX methods are strong-stability-preserving to maximize stable time-steps no matter whether the constraint is due to the implicitly integrated terms (diffusion) or the explicitly integrated ones (advection). From that point of view SSP IMEX methods with a non-trivial region of absolute monotonicity as defined by Theorem 1.1 are the most robust methods, because they allow achieving optimally large time-steps also at low resolution. We note here that while the region of absolute monotonicity has to be observed with respect to the explicitly integrated advection operator of the dynamical equations, stable time integration can be performed with step sizes falling outside it, if the restriction is due to the implicitly integrated diffusion terms. Thus, the class of optimal integrators for this kind of problem is probably larger than that of the SSP IMEX methods. However, none of the other time integration methods could significantly outperform them with respect to the time-steps achievable, and we have always found at least one case, where competing methods fell substantially short or even failed.
There is potential to further optimize the implementation of IMEX methods. The additional computational effort due to the implicit subpart is compensated for by accuracy and stability, but could be reduced in the future by replacing the solver for the linear equations associated with the arising generalized Poisson problem by a multigrid solver. In the Boussinesq approximation, additional solution of a Helmholtz equation is necessary instead. This widens the choice of fast solvers for the system of linear equations introduced through implicit time integration. The benefits expected from faster solvers would allow taking full advantage of the potential of method (33) that is implied by the large time steps reported in Table 8. Such an improvement would likewise be useful for the present problem to minimize the overhead by any of the implicit schemes in the regime where the time-step is limited by τ fluid rather than τ T . 1 + µ √ 2 − 2 µ + 3 µ 2 − 2 µ 2 √ 2 + 7/6 µ 3 √ 2 − 5/3 µ 3 π 2 1 − 2 µ + 2 µ 2 − 2/3 µ 3 π 1 − 4 µ + 8 µ 2 − 16/3 µ 3 Table 11 Values of g(θ, µ) for some θ, SSPRK(3,2) scheme (explicit scheme A from (33)), three point space discretization (8).  Table 12 Values of g(θ, µ) for some θ, SSPRK(3,2) (explicit scheme A from (33)), fourth order space discretization (9).
The point where the stability region intersects the negative real half-line is located at x ≈ −4.519.
The dissipativity analysis for A yields the amplification factors for the standard three-point space discretization (8) and the fourth order stencil (9). These amplification factors are evaluated at the points θ ∈ {0, π 4 , π 2 , π} in Tables 11  and 12, respectively.