The time step constraint in radiation hydrodynamics

TL;DR

This paper investigates the radiative time step constraint in explicit radiation hydrodynamics simulations, revealing it is governed by combined optically thick and thin contributions, and discusses potential mitigation strategies.

Contribution

It identifies the combined effect of optically thick and thin contributions on the radiative time step constraint and empirically determines their weighting in simulations, offering insights for improved modeling.

Findings

The radiative time step constraint is governed by the sum of optically thick and thin contributions.

In convective accretion disc simulations, the Deardorff term dominates enthalpy flux.

Mitigation strategies for the radiative time step problem are discussed.

Abstract

Explicit radiation hydrodynamic simulations of the atmospheres of massive stars and of convection in accretion discs around white dwarfs suffer from prohibitively short time steps due to radiation. This constraint is related to the cooling time rather than the radiative pressure, which also becomes important in hot stars and discs. We show that the radiative time step constraint is governed by the minimum of the sum of the optically thick and thin contributions rather than the smaller one of the two. In simulations with the Pencil Code, their weighting fractions are found empirically. In three-dimensional convective accretion disc simulations, the Deardorff term is found to be the main contributor to the enthalpy flux rather than the superadiabatic gradient. We conclude with a discussion of how the radiative time step problem could be mitigated in certain types of investigations.

Equations

$$\hat{\mbox{\boldmath$n$}}{}\bm{\cdot}\mbox{\boldmath$\nabla$}{}I=-\kappa\rho\,(I-S),$$

$${\textstyle{1\over 3}}(\ell\mbox{\boldmath$\nabla$}{})^{2}J=J-S,$$

$$\rho c_{\!p}{{\rm D}{}T\over{\rm D}{}t}-{{\rm D}{}p\over{\rm D}{}t}={\textstyle{4\pi\over 3}}\,\kappa\rho\,(\ell\mbox{\boldmath$\nabla$}{})^{2}J,$$

$$\lambda =\frac{c_{\gamma }}{\ell }\frac{k^2{\ell }^2/3}{1+k^2{\ell }^2/3}=\frac{c_{\gamma }{k}^2\ell /3}{1+k^2{\ell }^2/3}=\frac{\chi k^2}{1+k^2{\ell }^2/3},$$

$$c_{\gamma }=16{\sigma }_{\mathrm{SB}}{T}^3/\rho c_p$$

$$\lambda_{\rm exact}={c_{\gamma}\over\ell}\,\left(1-\frac{1}{k\ell}{\rm arccot}\,{}\frac{1}{k\ell}\right)={c_{\gamma}\over\ell}\,\left(1-\frac{\arctan k\ell}{k\ell}\right)\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil$\displaystyle#$\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil$\textstyle#$\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil$\scriptstyle#$\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil$\scriptscriptstyle#$\hfil\cr<\cr\sim\cr}}}}\lambda.$$

$${\tau }_{\mathrm{cool}}=\frac{1}{\chi k^2}+\frac{\ell }{c_{\gamma }}.$$

$$\delta t_{\mathrm{rad}}=C_{\mathrm{rad}}^{\mathrm{thick}}\frac{\delta z^2}{\chi D}+C_{\mathrm{rad}}^{\mathrm{thin}}\frac{\ell }{c_{\gamma }},$$

$$\delta t\le \min (\delta t_{\mathrm{rad}}),$$

$$\delta t_{\mathrm{rad}}^{\mathrm{Athena}}\propto \min (\text{mbox}\mathrm{Bo})\min (\delta z/c_{\mathrm{s}});$$

$$\delta t=\min \left(C_{\mathrm{rad}}^{\mathrm{thick}}\frac{\delta z^2}{\chi D}+C_{\mathrm{rad}}^{\mathrm{thin}}\frac{\ell }{c_{\gamma }},C_{\mathrm{CFL}}\frac{\delta z}{c_{\mathrm{s}}},C_{\mathrm{visc}}\frac{\delta z^2}{\nu D}\right).$$

$$\delta t\ne \min \left(C_{\mathrm{rad}}^{\mathrm{thick}}\frac{\delta z^2}{\chi D},C_{\mathrm{rad}}^{\mathrm{thin}}\frac{\ell }{c_{\gamma }},C_{\mathrm{CFL}}\frac{\delta z}{c_{\mathrm{s}}},C_{\mathrm{visc}}\frac{\delta z^2}{\nu D}\right).$$

$${{\rm D}{}\ln\rho\over{\rm D}{}t}=-\mbox{\boldmath$\nabla$}{}\bm{\cdot}\bm{u},$$

$$\rho{{\rm D}{}\bm{u}\over{\rm D}{}t}=-\mbox{\boldmath$\nabla$}{}p+\rho\mbox{\boldmath$g$}{}+{\rho\kappa\over c}\bm{F}_{\rm rad}+\mbox{\boldmath$\nabla$}{}\bm{\cdot}\mbox{\boldmath$\tau$}{},$$

$$\rho T{{\rm D}{}s\over{\rm D}{}t}={\cal H}-\mbox{\boldmath$\nabla$}{}\bm{\cdot}\bm{F}_{\rm rad}+\mbox{\boldmath$\tau$}{}\bm{:}\mbox{\boldmath$\nabla$}{}\bm{U},$$

$$\hat{\mbox{\boldmath$n$}}{}\bm{\cdot}\mbox{\boldmath$\nabla$}{}I=-\kappa\rho\,(I-S),\quad\bm{F}_{\rm rad}=\int_{4\pi}\hat{\mbox{\boldmath$n$}}{}I\,{\rm d}{}\Omega,\quad\mbox{\boldmath$\nabla$}{}\bm{\cdot}\bm{F}_{\rm rad}=\int_{4\pi}(I-S)\,{\rm d}{}\Omega,$$

$$\nu_{\rm shock}=C_{\rm shock}\delta x^{2}\langle-\mbox{\boldmath$\nabla$}{}\bm{\cdot}\bm{u}\rangle_{+},$$

$$\tau_{ij}=2\rho\nu{\sf S}_{ij}+\rho\nu_{\rm shock}\delta_{ij}\mbox{\boldmath$\nabla$}{}\bm{\cdot}\bm{u}.$$

$$T(z,t=0)=T_0+T_1\sin kz,$$

$$\kappa ={\kappa }_{\mathrm{es}}\approx 0.34{\mathrm{cm}}^2{\mathrm{g}}^{-1}.$$

$$g\to g_{\mathrm{eff}}=g-(\kappa /c)F_{\mathrm{rad}}=\mathrm{const}.$$

$$\frac{\mathrm{d}p}{\mathrm{d}z}=-\rho g_{\mathrm{eff}},\frac{\mathrm{d}T}{\mathrm{d}z}=-\frac{F_{\mathrm{rad}}}{K},$$

$$\frac{\mathrm{d}T}{\mathrm{d}p}=\frac{F_{\mathrm{rad}}}{K\rho g_{\mathrm{eff}}}.$$

$${\nabla }_{\mathrm{rad}}\equiv \frac{\mathrm{d}\ln T}{\mathrm{d}\ln p}=\frac{p{F}_{\mathrm{rad}}}{KT\rho g_{\mathrm{eff}}}=\frac{(c_p-c_v)F_{\mathrm{rad}}}{K{g}_{\mathrm{eff}}},$$

$$\nabla =\min ({\nabla }_{\mathrm{ad}},{\nabla }_{\mathrm{rad}}).$$

$$(\nabla-\nabla_{\rm rad})+\epsilon_{\ast}(\nabla-\nabla_{\rm rad})^{\xi}=0\qquad\mbox{(for $\nabla_{\rm rad}>\nabla_{\rm ad}$)}$$

$$\frac{1}{\kappa }=\frac{1}{{\kappa }_{{\mathrm{H}}^{-}}}+\frac{1}{{\kappa }_{\mathrm{Kr}}+{\kappa }_{\mathrm{es}}},$$

$${\kappa }_i={\kappa }_0(\rho /{\rho }_0{)}^{a_i}(T/T_0{)}^{b_i}$$

$$\mathcal{H}(z)={\left(\frac{3}{2}\Omega \right)}^2\frac{\dot{M}}{3\pi }\frac{\Theta (z_{\mathrm{heat}}-|z|)}{z_{\mathrm{heat}}},$$

$$\tau (z)={\int }_z^{\infty }\kappa (z^{\prime })\rho (z^{\prime })\mathrm{d}{z}^{\prime },$$

Full text

††articletype: GAFD Special issue on “Physics and Algorithms of the Pencil Code”

1 Introduction

Numerical simulations have long played an essential role in facilitating our understanding of hydrodynamic processes in astrophysics. The cost of such simulations is determined not only by the numerical mesh resolution, but also by the length of the time step. In hydrodynamics with explicit time-stepping, the maximum permissible time step decreases linearly with increasing spatial mesh resolution in such a way that the information that is passed from one time step to the next cannot propagate by more than roughly one mesh spacing $\delta x$. The length of the time step $\delta t$ is therefore of the order of $\delta x/c$, where $c$ is the speed of the fastest propagating mode (for example the speed of sound in subsonic compressible simulations). This is known as the Courant-Friedrichs-Lewy (CFL) condition (Courant et al., 1928) in ordinary hydrodynamics. This condition is known to change for diffusive processes with diffusivity $\chi$. In radiative flows, the time taken to propagate information from one mesh point to the next is often estimated based on the diffusion approximation, which gives \delta t\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\hfil\displaystyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\textstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptstyle#\hfil\cr<\cr\sim\cr}}}{\vbox{\offinterlineskip\halign{\hfil\scriptscriptstyle#\hfil\cr<\cr\sim\cr}}}}C_{\rm diff}\delta x^{2}/\chi (see, e.g., Castor, 2004), where $C_{\rm diff}$ is an empirical parameter, which Caunt and Korpi (2001) found to be $0.05$ in the context of magnetic diffusion. Thus, decreasing the mesh width by a factor of two implies a reduction of the time step by a factor of four. Indeed, Castor (2004) quotes this as the main reason against explicit radiation hydrodynamics in general. While this may be true for certain cases, one must recall that the diffusion approximation is valid only in the optically thick regime.

As one approaches the outer layers, the opacity decreases sharply and the mean free path becomes long compared to other typical scales in the system. If one were to continue using the optically thick approximation, the diffusivity would become excessively large. This would have severe consequences for the length of the time step. In reality, however, the optically thick approximation becomes invalid and no stringent time step constraint is expected to occur in the optically thin regime – at least not for solar-type stars, although this may change when more realistic opacities are invoked (Freytag et al., 2012). However, empirically we know that for hotter stars there can be layers in the proximity of the photosphere where the time step constraint can become rather stringent.

In the present paper, we will be concerned with radiation transport using what is nowadays often called long characteristics (Mihalas, 1978; Nordlund, 1982). Relativistic effects are ignored and the radiation field is propagated instantaneously across all rays without imposing any direct time step constraint. We refer to Finlator et al. (2009) for a treatment of time-dependent radiative transfer simulations of cosmological reionisation using long characteristics and to Pomraning (1979) for a discussion of the non-equilibrium Marshak wave problem. We also assume that the source function is just given by the Planck function and thus ignore the possibility that it depends on the mean intensity. This implies that scattering is treated as true absorption, as is commonly done (Freytag et al., 2012), but see the work of Skartlien (2000) for a detailed treatment of scattering. However, radiation interacts with the velocity field through radiation pressure and the temperature field through heating and cooling processes. While neither of these processes usually impose computationally prohibitive time step constraints in solar physics (Stein and Nordlund, 1989, 1998), a serious time step constraint (more stringent than the hydrodynamic time step) has been encountered empirically in numerical solutions of hot atmospheres where radiation pressure contributes to the hydrostatic equilibrium (Spiegel, 2006). One possibility is that the radiative pressure was responsible for the empirically determined short time step. More recently, the authors of the present work have encountered a similar time step constraint when solving the radiation hydrodynamics equations for hot accretion discs around white dwarfs. However, since the radiation pressure was not included in those solutions, it could not be held responsible here.

The fact that radiative time scales can be more restrictive than hydrodynamic ones, and would hence lead to a stricter time step constraint in explicit solvers, is not new (see, e.g., Castor, 2004; Davis et al., 2012). However, there seem to be conflicting statements regarding this problem. For example, Davis et al. (2012) derived a generalised CFL condition for the explicit radiation transfer solver in Athena (Stone et al., 2008) using short characteristics. However, instead of considering a general expression valid in both optically thick and thin regions, they approximated the radiative time step to be proportional to the time step due to the hydrodynamical CFL constraint locally in a grid zone (see section 2 for more details). This radiative time step is switched on whenever the optical depth per grid zone drops below unity anywhere in their computational domain. Such an assumption does not seem justified, as will be explained below. Furthermore, Davis et al. (2012) report on discrepancies and problems with lower optical depths, as well as very short time steps when the radiation energy density dominates the gas energy density. There is also a tendency to resort to semi-implicit (see, e.g., González et al., 2007, for the HERACLES code), or fully-implicit radiation transfer solvers (see, e.g., Stone and Norman 1992a, b; Stone et al. 1992 for the Zeus code; see also Jiang et al. 2012 for another module of Athena) to bypass the radiative time step problem. However, while implicit methods can avoid very small time steps, they are computationally more expensive. Freytag et al. (2012) proposed yet another time step constraint that is valid both in optically thick and thin regimes. However, no detailed study is presented. Their proposal qualitatively agrees with ours, but is different in that we suggest the additional presence of free parameters that have not been suggested or motivated before, and that alleviate the constraint in optically thin regions.

In view of the different and sometimes conflicting proposals for the radiative time step in explicit radiation hydrodynamics, there appears to be a need for a more a rigorous investigation. Our hope is that by diagnosing in more detail the radiative time step constraints in different situations, we would be in a better position to avoid or mitigate the problem of short time steps. One possibility might be to adopt certain changes in the model setup, while still being able to capture the essential physics. This will be discussed at the end of this paper.

2 Radiative cooling constraint

To quantify the expected time step constraint in radiation hydrodynamics, we begin by computing the cooling time. For that, we have to consider the radiation transport equation for the intensity I(\bm{x},t,\hat{\mbox{\boldmathn}}{}), where $\bm{x}$ is position, $t$ is time, and \hat{\mbox{\boldmathn}}{} is the direction of the ray. In the grey approximation, the radiation transport equation is

[TABLE]

where $\rho$ is the density, $\kappa$ is the opacity per unit mass, and $S(\bm{x},t)$ is the source function, which we will assume to be given by the Planck function, i.e., $S=(\sigma_{\rm SB}/\pi)T^{4}$, with $\sigma_{\rm SB}$ being the Stefan-Boltzmann constant, and $T$ is the temperature.

To gain insight into the nature of radiation in the optically thin and thick cases, it is useful to adopt a model where we can assume constant coefficients, which allows us to use Fourier transformation. We also adopt the Eddington approximation, where the moment expansion is closed by assuming the radiation pressure to be isotropic and given by ${\textstyle{1\over 3}}\delta_{ij}J$, where $J=\int_{4\pi}I\,{\rm d}{}\Omega/4\pi$ is the mean intensity and ${\rm d}{}\Omega$ is the differential over the solid angle. This yields (Edwards, 1990)

[TABLE]

where $\ell=(\kappa\rho)^{-1}$ is the photon mean-free path. Note that (\ell\mbox{\boldmath\nabla}{})^{2}=\ell^{2}\mbox{\boldmath\nabla}{}^{2} only when $\ell$ is constant in space. In the absence of any heating and cooling processes other than the negative radiative flux divergence, -{\mbox{\boldmath\nabla}{}\bm{\cdot}\bm{F}_{\rm rad}}, which is proportional to (\ell\mbox{\boldmath\nabla}{})^{2}J, the temperature evolution is governed by the radiative heat equation

[TABLE]

which is valid both in the optically thick and thin cases. Here, $p$ is the pressure, $c_{\!p}$ is the specific heat at constant pressure, and {\rm D}{}/{\rm D}{}t=\partial/\partial t+\bm{u}\bm{\cdot}\mbox{\boldmath\nabla}{} is the advective derivative. For the purpose of the present discussion, we assume $p={\rm const}{}$ and omit the ${\rm D}{}p/{\rm D}{}t$ term. Linearising equation (3) about a hydrostatic homogeneous equilibrium solution with $\bm{u}=\bm{0}$, $T={\rm const}{}$, and $\rho={\rm const}{}$, and assuming the solution to be proportional to $e^{{\rm i}\bm{k}\bm{\cdot}\bm{x}-\lambda t}$, where $\bm{k}$ is the wavevector, we find for the cooling or decay rate $\lambda$ the expression (Unno and Spiegel, 1966)

[TABLE]

where $k=|\bm{k}|$ is the wavenumber,

[TABLE]

is the characteristic velocity of photon diffusion (Barekat and Brandenburg, 2014), and $\chi=c_{\gamma}\ell/3$ is the radiative diffusivity. The quantity $c_{\gamma}$ is related to the radiative relaxation time $\ell/c_{\gamma}$ (equivalent to $q^{-1}$ of Unno and Spiegel 1966). It is smaller than the speed of light by roughly the ratio of radiative to thermal energies. Expression (4) has been obtained under the Eddington approximation and deviates only slightly from the exact expression obtained by Spiegel (1957), which can be written as

[TABLE]

The largest value of the ratio $\lambda/\lambda_{\rm exact}$ is 1.29 at $k\ell=2.53$.

The former expression (4) has the advantage that its inverse is the sum of two terms, allowing us to easily analyse the different regimes. When $\ell$ is small in the sense that $k\ell\ll 1$, the cooling rate is $\lambda\approx\chi k^{2}=c_{\gamma}\ell k^{2}/3$, which corresponds to the usual expression in the optically thick limit. On the other hand, in the optically thin limit, when $\ell$ is large ($k\ell\gg 1$), we have $\lambda\approx c_{\gamma}/\ell$, so cooling becomes independent of $k$. The resulting value of $\lambda$ is much smaller than the value of $c_{\gamma}\ell k^{2}/3$, which it would be if one continued using the expression for the optically thick case even though $k\ell\gg 1$. Thus, when using the correct expression, it may seem that radiation is less likely to be a limiting factor in the time step consideration. However, as we will show in this work, this is not necessarily true.

To arrive at an expression for the time step constraint due to radiation, we follow the idea of Freytag et al. (2012) and express the cooling time as the inverse of the cooling rate, $\tau_{\rm cool}=\lambda^{-1}$, where

[TABLE]

The first and second terms on the right-hand side of (7) characterise the contributions from the optically thick and thin parts, respectively, and $k$ should be replaced by the Nyquist wavenumber, $k_{\rm Ny}=\pi/\delta z$. Here and in the following sections, we restrict our attention to the vertical mesh spacing $\delta z$. However, in more than one dimension, the $k^{2}$ factor gains additional contributions $k_{x}^{2}=(\pi/\delta x)^{2}$ and $k_{y}^{2}=(\pi/\delta y)^{2}$ for the mesh spacings $\delta x$ and $\delta y$ in the $x$ and $y$ directions. Assuming $\delta x=\delta y=\delta z$, we would need to replace $k^{2}\to k_{x}^{2}+k_{y}^{2}+k_{z}^{2}=3(\pi/\delta z)^{2}$, or, more generally for $D$ dimensions, by $(\pi/\delta z)^{2}D$. In general, the direction with the finest mesh spacing will impose the strongest constraint.

The parts for the optically thick and thin regimes may contribute with different non-dimensional prefactors to the actual (empirically determined) time step constraint. Nevertheless, we will assume that the time step limit still depends on the sum of these two parts, just like in equation (7), because this ensures that in the optically thin regime, the much faster cooling rate relevant to the optically thick limit does not contribute. For a global simulation to be stable in both regimes, the time step can therefore not exceed the shortest value of the sum anywhere in the domain. Keeping this in mind, we define the radiative time step constraint to be

[TABLE]

such that the corresponding maximum permissible time step is given by

[TABLE]

where $C_{\rm rad}^{\rm thick}$ and $C_{\rm rad}^{\rm thin}$ are dimensionless coefficients. If we replace $k$ by the Nyquist wavenumber in $D$ dimensions, as done above, we have $C_{\rm rad}^{\rm thick}=1/\pi^{2}$, To estimate the value of $C_{\rm rad}^{\rm thin}$, we must consider the properties of the time stepping scheme. If we used a first order Euler scheme, then $\delta t\,c_{\gamma}/\ell<2$ is required for stability (Stoer and Bulirsch, 2002), i.e., $C_{\rm rad}^{\rm thin}<2$. Since this factor is larger than unity, it also has implications for the actual value of $C_{\rm rad}^{\rm thick}$, which would therefore be closer to $C_{\rm rad}^{\rm thick}=2/\pi^{2}\approx 0.2$. If we were to combine the Euler time stepping scheme with a less accurate second order discretisation of the Laplacian, the effective Nyquist wavenumber would only be $2/\delta z$ instead of $\pi/\delta z$, and one would obtain $C_{\rm rad}^{\rm thick}<2/2^{2}=0.5$. This is a well known result from the von Neumann stability analysis of the heat equation (Crank and Nicolson, 1947; Charney et al., 1950).

Note that $c_{\gamma}$ enters in both terms of equation (8). In the first term, it indeed plays the role of a photon diffusion velocity, but in the second term it represents a characteristic photon crossing velocity in an optically thin medium.

In direct contrast with the above calculation, we quote the radiative time step constraint used by Davis et al. (2012),

[TABLE]

see their equations (29) and (43), where $\mbox{\rm Bo}=16\,c_{\rm s}/c_{\gamma}$ is the local Boltzmann number and $c_{\rm s}$ is the sound speed. This time step constraint is used in the explicit radiation transfer module of Athena. The time step is here explicitly proportional to the usual CFL condition. Thus, equation (10) does not properly account for the transition from the optically thick to thin regions, as mentioned in the introduction. Note that the Boltzmann number is commonly invoked in radiative flows (see, e.g., Castor, 2004; Davis et al., 2012), and it is a measure of when radiation becomes important in a problem.

Equation (7) was also the basis of the time step constraint used in the CO5BOLD code (Freytag et al., 2012), which is used to model solar and stellar surface convection. In that code, radiation transport is also treated explicitly (although there are modules allowing for semi-implicit or fully implicit treatments). According to the CO5BOLD user manual111See http://www.astro.uu.se/~bf/co5bold_main.html., there is only one prefactor, which they refer to as the “radiative Courant” factor and they recommend it to be adjusted by trial and error. In this work, we introduce two separate coefficients $C_{\rm rad}^{\rm thick}$ and $C_{\rm rad}^{\rm thin}$, which can be thought of as the radiative Courant-like factors in the optically thick and thin regions, respectively. Both coefficients depend on the time stepping scheme, but the first one also depends on the spatial discretisation scheme.

It is useful to reflect again on the somewhat unusual form of equation (9) as a time step constraint, because the usual CFL condition is formulated in terms of the shortest one of several constraints, e.g., $\delta t=\min(\delta z/c_{\rm s},\,\delta z^{2}/\chi D)$. Alternatively, one could express the inverse time step as the sum of the inverse of the contributions. In (8), by contrast, $\delta t_{\rm rad}$ itself is determined as a sum of two time steps. In full radiation hydrodynamics, the maximum permissible time step will therefore be the minimum of $\delta t_{\rm rad}$, as given by equation (8), and the usual CFL and viscous constraints, i.e.,

[TABLE]

For the Pencil Code with its default third order time stepping scheme, $C_{\rm CFL}=0.9$ is the usual CFL number, $C_{\rm visc}=0.25$ is the viscous time step constraint, with $\nu$ being the kinematic viscosity and $D$ the dimensionality. In our present one-dimensional calculation, we have $D=1$. In three dimensions, $C_{\rm visc}/D=0.08$ is in fact slightly larger than the value $C_{\rm diff}=0.05$ quoted by Caunt and Korpi (2001). Again, the dimension $D$ enters because the discretised form of the second derivative, relevant in the optically thick formulation, has a larger coefficient at the center point where the derivative is evaluated.222For a second order discretisation, for example, we have $\nabla^{2}f_{i}=(f_{i+1}-2f_{i}+f_{i-1})/\delta x^{2}$ in one dimension, but $\nabla^{2}f_{ij}=(f_{i+1\,j}+f_{i\,j+1}-4f_{ij}+f_{i-1\,j}+f_{i\,j-1})/\delta x^{2}$ in two dimensions, so at the center point, $f_{ij}$, the coefficient increases from $2/\delta x^{2}$ to $4/\delta x^{2}$ in two dimensions, and to $6/\delta x^{2}$ in three dimensions. This applies analogously also to the sixth order discretisation used in the Pencil Code, except that the coefficient now increases to $49/18\approx 2.72$ instead of 2 per direction as for the second order case. Also, for a given function $f$, the value of $-(\nabla^{2}f)/f$ depends on the numerical scheme. For a checkerboard pattern of $f$ (e.g., an alternating sequence with $-1$, $+1$, etc, in one spatial dimension), using a second order scheme, the value is $-4/\delta z^{2}$ per direction, while the analytic value is $-\pi^{2}/\delta z^{2}\equiv k_{\rm Ny}^{2}$, which is more than twice as much.

We emphasise that there is no analogy in how the optically thick and thin contributions enter into the time step constraint and how the usual CFL and viscous constraints enter. This becomes strikingly clear by stating

[TABLE]

Later in this paper, we will see examples where either $\delta z^{2}/\chi$ or $\ell/c_{\gamma}$ may be very small, and yet, neither of those affect the time step if the other term is large.

The goal of the present paper is to test the validity of (8) in the case when radiation is treated with long characteristics. We also compare with the usual CFL condition, where the time step is constrained by $\delta t_{\rm s}=C_{\rm CFL}\delta z/c_{\rm s}$. Since the role of $c_{\gamma}$ is not entirely clear, especially when radiation pressure also enters the problem, we ask whether a similarly defined quantity $\delta t_{\gamma}=C_{\gamma}\delta z/c_{\gamma}$, as suggested by the work of Davis et al. (2012), might constrain the time step further, even though it does not explicitly feature in equation (11). Our numerical experiments reported below show that $\delta t_{\gamma}$ itself does not constrain the time step, although the ratio $c_{\rm s}/c_{\gamma}$ may still play an important role; see discussion in section 5. This can be better understood by recalling that the Boltzmann number Bo is proportional to $c_{\rm s}/c_{\gamma}$. Thus, a smaller value of Bo signifies that the energy transport is radiation dominated and vice versa.

If $c_{\gamma}$ were to enter the time step constraint directly through a quantity $t_{\gamma}$, it would be natural to expect that $C_{\gamma}$ would be of the order of $C_{\rm CFL}$. Hence, we use $C_{\gamma}=C_{\rm CFL}=0.9$ in our plots below. Regarding the value of $C_{\rm rad}^{\rm thick}$, we expect it to be comparable to $C_{\rm visc}=0.25$, but our experiments reported below seem to be consistent with a slightly smaller value of $C_{\rm rad}^{\rm thick}=0.2$, so we will use that value in all the corresponding plots shown below. Finally, regarding the value of $C_{\rm rad}^{\rm thin}$, it is important to note that it enters without a $\delta z$ term; see second term in (8). In order to obtain a preliminary estimate, we empirically test the radiative time step constraint in the optically thin case by using a one-dimensional model. In this way, we find that $C_{\rm rad}^{\rm thin}\approx 4$ and we use this value in the plots shown below; see section 3.3 for details. We find that the values of the various coefficients quoted here are consistent with the empirically determined ones, as discussed in detail in section 4.3 and table LABEL:Tsum.

In this paper, we discuss two distinct models where severe time step constraints have been encountered. One is the model of Spiegel (2006) and the other is a local model of an accretion disc, similar to that of Coleman et al. (2018). Radiation pressure is included in the former, but not in the latter. It will turn out that the time step constraints are quite different from each other in the two cases, although this difference is not explicitly linked to the presence or absence of radiation pressure. Nevertheless, when the radiation pressure becomes extremely large, it can in principal also restrict the time step. We discuss this possibility at the end of our penultimate section 5.

Finally, let us note that the radiation magnetohydrodynamic shearing-box simulations of accretion discs by Coleman et al. (2018) were carried out using the Zeus code. Time step problems were probably not encountered in this case, as Zeus employs an implicit radiation transfer solver.

3 Our models

3.1 The basic equations

We consider here the nonrelativistic radiation hydrodynamics equations solved by default in the Pencil Code. Note, however, that we work in what is referred to as the static diffusion limit. Defining $\beta=u/c$, where $u$ is the characteristic velocity of the system and $c$ the speed of light, the static diffusion limit is valid when $\beta\tau\ll 1$ in optically thick regions where $\tau\equiv\int\rho\kappa\,{\rm d}{}z\gg 1$; see Krumholz et al. (2007) for a detailed discussion of the various regimes in radiation hydrodynamics.

The basic dependent variables are the logarithmic density $\ln\rho$, the velocity $\bm{u}$, and the specific entropy $s$, which obey the equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where \mbox{\boldmathg}{}=(0,0,-g) is the gravitational acceleration in Cartesian coordinates $(x,y,z)$, $\bm{F}_{\rm rad}$ is the radiative flux, \mbox{\boldmath\tau}{}=2\rho\nu\mbox{\boldmath{\sf S}}{} is the stress tensor if there is just shear viscosity, {\sf S}_{ij}={\textstyle{1\over 2}}(\partial_{i}u_{j}+\partial_{j}u_{i})-{\textstyle{1\over 3}}\delta_{ij}\mbox{\boldmath\nabla}{}\bm{\cdot}\bm{u} are the components of the traceless rate-of-strain tensor, ${\cal H}$ is a heating function to be specified below, \mbox{\boldmath\tau}{}\bm{:}\mbox{\boldmath\nabla}{}\bm{U}\equiv\tau_{ij}\partial u_{i}/\partial x_{j} is the viscous heating term, and equation (1) has been restated in (16) for completeness. The specific entropy is related to pressure and density via ${\rm D}{}s=c_{\!v}{\rm D}{}\ln p-c_{\!p}{\rm D}{}\ln\rho$, where $c_{\!v}$ is the specific heat at constant volume. In the model for the upper layers of a star, we have $g={\rm const}{}$, while in our accretion disc model with local Keplerian angular velocity $\varOmega$ (not to be confused with the solid angle $\Omega$), $g=\varOmega^{2}z$ increases linearly with the distance $z$ from the midplane. Its sign is such that gravity pulls toward the midplane from above and below $z=0$.

The equations of radiation transport (16) have been implemented into the Pencil Code by Heinemann et al. (2006). This implementation was then used in the solar context (Heinemann et al., 2007) and later also in more idealised problems (Barekat and Brandenburg, 2014). In the outer layers of the Sun, partial ionisation is also important, so one needs to solve the Saha equation, for which temperature needs to be known. It is then advantageous to use $\ln T$ as the dependent variable instead of $s$. Again, this implementation into the Pencil Code goes back to the work of (Heinemann et al., 2007), and more idealised models with ionisation and radiation have been considered by Bhat and Brandenburg (2016).

In the presence of shocks, it is often useful to increase the viscosity locally in those regions where the velocity converges, i.e., where the flow divergence is negative or \mbox{\boldmath\nabla}{}\bm{\cdot}\bm{u}<0. This approach goes back to von Neumann and Richtmyer (1950). In practice, one defines the shock viscosity as

[TABLE]

where $\langle...\rangle_{+}$ denotes a running five point average over all positive arguments, $C_{\rm shock}$ is the dimensionless coefficient of the shock viscosity, and $\delta x=\delta y=\delta z$ are the mesh spacings in all three directions. This shock viscosity is applied as a bulk viscosity, i.e., $\tau$ in equations (14) and (15) is given by

[TABLE]

At the end of this paper, we present one case where a shock viscosity is applied and show how it affects the time step (figures 10 and 11 of section 4.4).

3.2 Numerical treatment within the Pencil Code

In the Pencil Code, all derivatives are usually approximated by sixth order finite differences; see Brandenburg (2003) for details. The third-order time stepping scheme of Williamson (1980) is used. Thus, during each of the three substeps, the right-hand side of the equations is evaluated three times. The code is not explicitly conservative, but mass, energy, and momentum are conserved within the discretisation error of the scheme, so we can use their conservation properties to gauge the accuracy of the scheme. Normally, when doubling the resolution, the error decreases by a factor of $2^{6}=64$, as was demonstrated by Babkovskaia et al. (2011).

The radiation transport is solved using long characteristics. A detailed account of its implementation has been given by Heinemann et al. (2006). For a brief description of the most important numerical aspects, we refer the reader to section 2.4 of Barekat and Brandenburg (2014). Instead of integrating along a geometric line segment ${\rm d}{}l$, we integrate (1) over optical depth ${\rm d}{}\tau=-\kappa\rho\,{\rm d}{}l$ for a number of rays with direction \hat{\mbox{\boldmathn}}{}, such that ${\rm d}{}I/{\rm d}{}\tau=I-S$. In formulating the boundary conditions, we distinguish between two types of rays. For rays that are perfectly horizontal, we assume periodicity. For all other rays, we normally assume that no radiation enters at the outer boundary of the domain (Mihalas, 1978). However, to reproduce the analytic solution for an infinitely extended domain, we can take the layers beyond the computational domain into account if the flux is assumed to be known and the solution is in radiative equilibrium. This is explained in appendix A.

The code is parallelised by splitting the problem into two local ones and a nonlocal one in between. The latter requires interprocessor communication. The local problems are computationally intensive, while the nonlocal problem does not involve any computations and is therefore quite fast. On each processor, one first solves (1) along each ray to compute the intrinsic intensity increment within each processor as a function of optical depth. In the second step, the increments of intrinsic intensity and optical depth are communicated to the neighbouring processors. In the third and final step, these increments are used to construct the total intensity within each processor. For most of the calculations presented here, we use just the vertical or $z$ direction, which corresponds to two ray directions, namely for upward and downward propagating radiation.

The radiative time step problem was never addressed in this code, because in the applications of Heinemann et al. (2007) to solar convection and sunspot formation, no severe time step constraints were encountered. Doing this quantitatively is the main purpose of the present paper.

3.3 Radiative time steps from an unstratified model

To have a preliminary idea about the values of $C_{\rm rad}^{\rm thick}$ and $C_{\rm rad}^{\rm thin}$, we perform simple one-dimensional experiments by solving equations (1) and (3) for constant $\rho$ and $\kappa$ (here with $\ell\equiv(\kappa\rho)^{-1}=1\,{\rm Mm}$) for an unstratified model. We adopt periodic boundary conditions for $I$ and $T$ and look for the decay of a sinusoidal temperature perturbation of the form

[TABLE]

where $T_{0}=10^{6}\,{\rm K}$ is chosen hot enough so that radiation provides the most restrictive time step constraint and $T_{1}=10^{5}\,{\rm K}$. The density is chosen to be $10^{-9}\,{\rm g}\,{\rm cm}^{-3}$. The value of $c_{\gamma}$ is then $2.2\times 10^{10}\,{\rm km}\,{\rm s}^{-1}$. This exceeds the speed of light, but that does not matter in this non-relativistic computation, and it helps making sure the acoustic time step remains unimportant under all circumstances. Thus, in the optically thin case, the cooling time $\ell/c_{\gamma}$ is much smaller than the sound wave crossing time, $\delta z/c_{\rm s}$.

Here we decided to retain equations (13) and (14), which implies that sound waves can equilibrate pressure perturbations on an acoustic time scale. This has another interesting side effect in that the temperature perturbations then tend to cause corresponding density variations rather than pressure variations. This justifies in hindsight that the ${\rm D}{}p/{\rm D}{}t$ term in (3) could then be omitted. This, in turn, justifies the presence of $c_{\!p}$ instead of $c_{\!v}$ in the definition of $c_{\gamma}$ in (5); see appendix B for details.

We vary $k$ in the range from $10^{-4}$ to $10\,{\rm Mm}^{-1}$ and choose $N=128$ mesh points, so the domain mesh widths are $\delta z=2\pi/kN$, which varies between approximately $500$ and $0.005\,{\rm Mm}$, respectively. In all cases, we determine the maximum permissible time step $\delta t$. An example of the time evolution of the maximum temperature $T_{\max}$ versus time is presented in figure 1, which show that $T_{\max}$ develops runaway when the time step is too long. In this case, we find that $C^{\rm thin}_{\rm rad}$ or $\delta t\,c_{\gamma}/\ell\approx 4.3300$ is the approximate borderline. Already 4.3302 is too large.

In figure 2, we show the resulting maximum permissible time step, $\delta t$, in two different normalisations related to the two terms of equation (7). When $\delta t$ is normalised by $\ell/c_{\gamma}$, which is a constant in this model, we see that $\delta t$ increases quadratically with $\delta z$ for large values, but is independent of $\delta z$ for small values – saturating at a value of about $3$. This is in perfect agreement with equation (8). On the other hand, when $\delta t$ is normalised by $\delta z^{2}/\chi$, it levels off at a value of around 0.3. Normalising instead by $(\chi k_{\rm Ny}^{2})^{-1}$, we see that the maximally permissible time step levels off at a value of about $3$, i.e., same as that for the normalisation by $\ell/c_{\gamma}$. This may justify the use of a single coefficient for the radiation time step, as done by Freytag et al. (2012). However, since the values of the two coefficients depend on the numerical scheme, we retain the two independently.

For small $\delta z$ (optically thin limit), we see that $\delta t\,c_{\gamma}/\ell\approx 3$–$4$, while for large $\delta z$ (optically thick limit), we obtain $\delta t\,\chi/\delta z^{2}\approx 0.3$. These, then, would be the recommended values of $C_{\rm rad}^{\rm thin}$ and $C_{\rm rad}^{\rm thick}$, respectively. However, throughout this work, we adopt a more conservative value for the latter, $C_{\rm rad}^{\rm thick}=0.2$, which is later found to be necessary for numerical stability in some stratified cases, while $C_{\rm rad}^{\rm thin}\approx 4$ turns out to work well in the cases presented below.

Let us now comment on the spatial properties of the solution in the optically thick and thin cases when the numerical instability develops, i.e., when the time step is too long. For $\delta z/\ell\gg 1$, i.e., when the time step is constrained by $\delta z^{2}/\chi$, the numerical instability develops its fastest growing mode uniformly over the domain at the mesh scale, i.e., at $k=k_{\rm Ny}$. This does indeed correspond to a checkerboard pattern in two dimensions, as was assumed in our analysis above. For $\delta z/\ell\ll 1$, on the other hand, the fastest growing mode also develops at $k=k_{\rm Ny}$, but nonuniformly and preferentially in those locations that are hotter.

We have stated in section 2 that the time step in the optically thick regime scales with the dimension $D$. This is demonstrated in table LABEL:Tdir, where we compare the maximum permissible $\delta t$ in units of $\delta z^{2}/\chi$ in multidimensional domains ($D=1$, $2$, and $3$) for different numbers of rays, $n_{\rm ray}$. In those cases, the mesh spacing is the same in all directions. We see that the time step is approximately inversely proportional to $D$ if we restrict ourselves to rays along the coordinate direction, i.e., $n_{\rm ray}=2$, $4$, or $6$. In that case, $\delta t\,\chi/\delta z^{2}\approx 0.38/D$. Interestingly, when more directions are included ($n_{\rm ray}=14$ or $22$), the minimal timestep becomes longer again. This is because the radiative flux divergence is calculated as an angular integral over all directions. However, the diagonal directions do not contribute for a checkerboard pattern. Therefore, the radiative flux divergence decreases for larger $n_{\rm ray}$, which alleviates the time step constraint correspondingly. It turns out that the extra cost associated with the use of more rays is easily being outweighted by being able to use a longer time step. The increased accuracy obtained by using more rays comes therefore effectively at no extra cost.

3.4 Stratification of hot stellar surface layers

Later in this paper, we address the question which time step constraint plays a role in which layers of certain stars. For that purpose, we need a simple model for stellar surface layers. We therefore present here the relevant equations that can simply be solved by integration. The same solutions can also be obtained using explicit time integration with the Pencil Code.

In hot stellar surface layers, the radiation pressure plays an important role. This is typically also the regime in which the electron scattering opacity is important (Frank et al., 1992), i.e.,

[TABLE]

Since this is a constant, and since the radiative flux $F_{\rm rad}$ is also constant in radiative equilibrium, we just have to replace gravity by the effective one, which is then still a positive constant, i.e.,

[TABLE]

To obtain a hydrostatic reference solution, we integrate the equations of hydrostatic and radiative equilibrium,

[TABLE]

with $F_{\rm rad}=\sigma_{\rm SB}T_{\rm eff}^{4}={\rm const}{}$, where $T_{\rm eff}$ is the effective temperature of the stellar surface and $K=\rho c_{\!p}\chi$. Note that the equation for ${\rm d}{}T/{\rm d}{}z$ applies even in the optically thin case provided the system is in thermal equilibrium. This is because in equilibrium we have \mbox{\boldmath\nabla}{}\bm{\cdot}\bm{F}_{\rm rad}=0, and therefore $J=S$, just like in the optically thick case. In the time-dependent case, however, as discussed above, the optically thick and thin cases are quite different from each other. It is only in that case that a time step constraint has to be obeyed.

Next, it is convenient to divide the two equations in (22) by each other, so we obtain,

[TABLE]

In astrophysics, the symbol $\nabla$ is used to denote the double-logarithmic temperature gradient, so $\nabla={\rm d}{}\ln T/{\rm d}{}\ln p$, and for the present radiative equilibrium solution, the radiative temperature gradient is denoted by $\nabla_{\rm rad}$. Thus, we have

[TABLE]

where we have used $p=(c_{\!p}-c_{\!v})\,T\rho$ for the equation of state of a perfect gas and $c_{\!v}$ is the specific heat at constant volume.

The vertical specific entropy gradient is ${\rm d}{}(s/c_{\!p})/{\rm d}{}\ln p=\nabla-\nabla_{\rm ad}$, so the Schwarzschild criterion for convective stability (positive outward gradient of $s$) is $\nabla<\nabla_{\rm ad}$. Convection occurs when $\nabla>\nabla_{\rm ad}$. This mixes the fluid, so $s$ becomes uniform and one must replace $\nabla$ by $\nabla_{\rm ad}=1-1/\gamma=0.4$, which is the value for a monatomic gas with $\gamma=c_{\!p}/c_{\!v}=5/3$. Thus, the local double-logarithmic temperature gradient can be written as (see, e.g., Kähler, 1972)

[TABLE]

In more realistic mixing length descriptions of stellar convection, this relation is to be replaced by a smooth transition between the two states; see Vitense (1953) for the original formulation, which corresponds to finding a solution to the equation

[TABLE]

with $\xi=3/2$, $\epsilon_{\ast}=\iota c_{\rm s}^{3}/\chi g$, and $\iota$ being a coefficient of the order of unity; see figure 3 of Brandenburg (2016) for a comparison of solutions for $\xi=3/2$ and $\xi=1$.

For lower temperatures, the opacity given by (20) is no longer valid. A more general representation is given in terms of combinations of Kramers-type opacities

[TABLE]

where we use

[TABLE]

with $i={\rm Kr}$ or ${\rm H^{-}}$. For our hot stellar surface models, we use $\kappa_{0}=10^{4}\,{\rm cm}^{2}\,{\rm g}^{-1}$, $\rho_{0}=10^{-5}\,{\rm g}\,{\rm cm}^{-3}$, $T_{0}=13,000\,{\rm K}$, and $a_{\rm Kr}=1$, $b_{\rm Kr}=-3.5$ for the Kramers opacity $\kappa_{\rm Kr}$, relevant for the deeper layers in the star, and $a_{\rm H^{-}}=0.5$, $b_{\rm H^{-}}=18$ for the H- opacity $\kappa_{\rm H^{-}}$ in the layers just beneath the photosphere. These were also the coefficients used by Brandenburg (2016). The same opacity prescription will also be used in our accretion disc models described below, but with $\kappa_{0}=2\times 10^{4}\,{\rm cm}^{2}\,{\rm g}^{-1}$, $T_{0}=20,000\,{\rm K}$, $a_{\rm H^{-}}=1$ and $b_{\rm H^{-}}=4$.

To construct a solution, we assume $T_{\rm eff}$ and $g_{\rm eff}$ as given. For grey atmospheres, the temperature of the atmosphere far above the photosphere is $T_{0}=T_{\rm eff}/2^{1/4}\approx 0.84\,T_{\rm eff}$; see Stix (2002) for a text book. Thus, we integrate from the top downward using $\ln p$ as the independent variable starting with a sufficiently low value. For single power law opacities, such as a simple Kramers-type opacity with fixed exponents $a$ and $b$, this integration can be done analytically (see appendix A of Brandenburg, 2016), but here we use more complicated opacities and do the integration numerically. At each height, we solve equation (26) for $\nabla$ to determine the temperature gradient when $\nabla_{\rm rad}>\nabla_{\rm ad}$, in which case we replace $\nabla$ by $\nabla_{\rm ad}$.

3.5 Accretion disc model

Accretion discs are generally heated by turbulent dissipation. To model this realistically, we would need to simulate accretion disc turbulence through a combination of (i) the magneto-rotational instability (MRI) to generate turbulence from a given magnetic field (Balbus and Hawley, 1991) and (ii) the dynamo instability to regenerate the required magnetic field (Brandenburg et al., 1995). To simplify matters, and to avoid modelling the MRI and the dynamo, we assume instead a prescribed heating function as a function of $z$,

[TABLE]

where we have used the usual parameterisation for accretion discs in terms of mass accretion rate $\dot{M}$ and local Keplerian angular velocity $\varOmega$ (Frank et al., 1992). In addition, we have assumed a vertical profile $\varTheta(z_{\rm heat}-|z|)$, which is unity inside the disc, $|z|\leq z_{\rm heat}$, and zero outside. In practice, we initiate our simulations by using an isothermal hydrostatic stratification, where $g=\varOmega^{2}z$ is the vertical gravity and $\ln(\rho/\rho_{0})=-z^{2}/2H_{\!p}^{2}$, with $H_{\!p}=c_{\rm s}/\varOmega$ being the pressure scale height, $c_{\rm s}=\sqrt{{\cal R}T/\mu}$ the isothermal sound speed, ${\cal R}$ the universal gas constant, $\rho_{0}$ the density at $z=0$, and $\mu$ the mean atomic weight. Radiation then causes the outer layers to cool until thermal equilibrium is achieved. Alternatively, one could start with a thermal hydrostatic equilibrium that is computed using equations (22)–(24), except that now $F_{\rm rad}\neq{\rm const}{}$. Hence, this has to be obtained by integrating ${\rm d}{}F_{\rm rad}/{\rm d}{}z={\cal H}(z)$ with $F_{\rm rad}(z=0)=0$ as a boundary condition at the midplane.

For the purpose of understanding the radiative time step constraint, we will not be concerned with convection in these simulations with the Pencil Code. Instead, we use a one-dimensional model in which the flow is either up or down, but not both, so no return flow and no convection are possible. Convection will, however, be discussed briefly in section 4.4, where we illuminate in more detail the properties of discs around white dwarfs. In such discs, the total flux is given by $F_{\rm tot}=F_{\rm rad}+F_{\rm conv}$, and it is found that a significant fraction of the convective flux $F_{\rm conv}$ is independent of the superadiabatic gradient, as is normally assumed in standard mixing length theory. This is particularly interesting in view of recent suggestions (Brandenburg, 2016) that even in solar and stellar convection, the convective flux may have a significant contribution from what is known as the Deardorff (1966, 1972) term. Numerical evidence for this flux was found in simulations of stellar convection (Käpylä et al., 2017, 2019), and the convective simulations of accretion discs in this work will provide further evidence for this.

Note that we have not made any attempt to smoothen the abrupt change in the heating function of (29) and have rather regarded this feature as an advantage. This is because it allows us to see whether this profile results in a similarly abrupt transition in the resulting vertical profiles of temperature or radiative flux. As will be demonstrated in section 4.3, the temperature profile turns out to be smooth, suggesting that no artefacts of the heating profile are introduced into the model. The total flux $F_{\rm tot}$, on the other hand, shows a sharp first derivative in the steady state. Again, this is not a problem, since it serves as a convenient “marker” of where the heating stops in the vertical profiles of various fluxes. We will return to this in section 4.4.

Both in the stellar and in the accretion disc models, we determine the location of the photosphere as the point where the optical depth,

[TABLE]

is unity. Here, $z\to\infty$ corresponds to a location far away from the disc or the star, but for the lower disc plane, the integral would need to go from $-\infty$ to $z$ instead.

4 Results

4.1 A hot stellar surface layer

We adopt here one of the stellar surface layer models of the unpublished work of Brandenburg and Spiegel, who considered a star of solar mass $M=M_{\odot}=2\times 10^{33}\,{\rm g}$, solar radius $R=R_{\odot}=7\times 10^{10}\,{\rm cm}$ and $g=g_{\odot}=GM_{\odot}/R_{\odot}^{2}=2.7\times 10^{4}\,{\rm cm}\,{\rm s}^{-2}$ (the solar value), but with a luminosity $L$ that is $2\times 10^{4}$ times the solar value $L_{\odot}=4\times 10^{33}\,{\rm erg}\,{\rm s}^{-1}$. The radiative flux is given by $F_{\rm rad}=L/(4\pi R^{2})=1.3\times 10^{15}\,{\rm erg}\,{\rm cm}^{-2}\,{\rm s}^{-1}$ and the effective temperature $T_{\rm eff}=69,000\,{\rm K}$. We solve the time-dependent equations (13)–(16) for the surface layers of the above star using the Pencil Code, with initial $\bm{u}=0$. We also put ${\cal H}=0$, $\nu=10^{12}\,{\rm cm}^{2}\,{\rm s}^{-1}$, and $C_{\rm shock}=3$. We refer to appendix A regarding the boundary conditions imposed on the intensity, and appendix C for those on the hydrodynamic variables. Our model has a depth of $60\,{\rm Mm}$ and uses 256 uniformly spaced mesh points. The $\tau=1$ surface is roughly in the middle of the domain, which we define to be at $z=0$. The temperature then varies between $110,000\,{\rm K}$ at $z=-30\,{\rm Mm}$ and $62,000\,{\rm K}$ at $z=30\,{\rm Mm}$. The sound speed varies between $c_{\rm s}=55\,{\rm km}\,{\rm s}^{-1}$ at the bottom and $c_{\rm s}=42\,{\rm km}\,{\rm s}^{-1}$ at the top, while $c_{\gamma}=55\,{\rm km}\,{\rm s}^{-1}$ at the bottom and about $10^{8}\,{\rm km}\,{\rm s}^{-1}$ at the top. The domain has a density contrast of $\rho_{\max}/\rho_{\min}\approx 300$, so the number of density scale heights is $\Delta\ln\rho=\ln 300\approx 6$. Note that, for this case, we have used $\kappa=\kappa_{\rm es}={\rm const}{}$ (i.e., in figure 3 only) instead of expression (27), but the difference would be minor.

In figure 3, we show the steady state temperature stratification for the above star (left panel), and the various candidates that could play a role in constraining the time step (right panel). The total duration of this simulation was $56.4\,{\rm s}$. Note that this figure, as well as the rest of the figures obtained from Pencil Code simulations in this paper, have been plotted at the final time of the corresponding simulation (by which the solution has attained a steady state), unless otherwise mentioned. We also compare with the empirically determined maximum permissible time step $\delta t$ obtained in the Pencil Code simulations, which turns out to be $4\,{\rm ms}$ for this case. Note that every simulation is run with a constant, pre-assigned time step. We recall that, in any time-dependent calculation, regardless of whether or not a steady state solution is reached, the simulation may become numerically unstable if the time step is too long. One can then keep decreasing the supplied time step (to the desired accuracy) in order to determine the limiting time step for that particular problem, namely $\delta t$ in our notation.

Interestingly, we find in figure 3 that $\delta t$ is very close to the line for $\delta t_{\rm rad}^{\rm thin}=4\,\ell/c_{\gamma}$. The line for $\delta t_{\rm rad}^{\rm thick}=0.2\,\delta z^{2}/\chi$ is well below the actual $\delta t$, which is consistent with (11) in that a term proportional to $\delta t_{\rm rad}^{\rm thick}$ does not constrain $\delta t$ in isolation. Similarly, $\delta t_{\gamma}=0.9\,\delta z/c_{\gamma}$ is not motivated to be a possible candidate for constraining the time step in this case. The acoustic time step $\delta t_{\rm s}=0.9\,\delta z/c_{\rm s}$ is much larger than the radiative one and therefore unimportant in this case. It is important to note here that the right-hand side of (12) would have predicted $\delta t_{\rm rad}^{\rm thick}$ to be the relevant time step in the problem. However, in reality it turns out to be $\delta t_{\rm rad}^{\rm thin}$, as also supported by our numerical simulations.

Spiegel and Tao (1999) motivated the interest in studying the hydrodynamics of hot stars by referring to O. Struve for having discovered large line widths, which, in turn, could hint at the existence of turbulence in the atmospheres of those stars. They associated this line broadening with photofluid instabilities, which are possible even when the radiative acceleration is still below the gravitational one. In the present case, the radiative acceleration is $0.53$ times the gravitational one, and may even exceed it at higher temperatures. This could lead to photoconvection (Prendergast and Spiegel, 1973; Spiegel, 1977; Spiegel and Tao, 1999), which Brandenburg and Spiegel attempted to study with the Pencil Code, but suspended this project because of prohibitively short time steps. In recent years, similar studies have been performed in the context of radiation-driven stellar winds (Owocki and Sundqvist, 2018; Sundqvist et al., 2018). No time step problems, have been reported in these studies. However, we point out that the radiation is treated differently there, as they consider two-dimensional radiation line-transport. Also, the time step is chosen such that it is the minimum of a fixed time and a variable 1/3 of the Courant time.

4.2 Expected time step constraints for stellar surface layers

We now consider solutions of the time-independent equations (21)–(27) for stellar surface layers (see section 3.4 for the method). So no time step constraint applies, but we can still use this model to predict what the time step constraint would be in a time-dependent simulation. In figure 4, we plot $T(\ln p)$ and $s(\ln p)$ for five values of $T_{\rm eff}=5000\,{\rm K},7000\,{\rm K},10000\,{\rm K},20000\,{\rm K}$ and $70000\,{\rm K}$. We use $g=g_{\odot}$ for the following plots. In figure 5, we plot the two contributions to the radiative time step, as well as the acoustic one,

[TABLE]

together with the total radiative one $\delta t_{\rm rad}=\delta t_{\rm rad}^{\rm thick}+\delta t_{\rm rad}^{\rm thin}$ for $\delta z=0.05H_{\!p}$, where $H_{\!p}=(c_{\rm p}-c_{\rm v})T/g$ is the local pressure scale-height in the star. In all of our models, the temperature above the photosphere reaches a constant; see figure 4(a). This is because the physics of realistic coronal heating and cooling is not included in our simple model; see the papers by Bingert and Peter (2011) and Bourdin et al. (2013) for realistic coronal modelling with the Pencil Code, using a setup originally developed by Gudiksen and Nordlund (2002, 2005a, 2005b).

Next, we discuss the constraints on the resulting time step, that would limit a time-dependent calculation of the same problem. In figure 5, at small values of $\ln p$, corresponding to locations above the photosphere, we have $\delta t_{\rm rad}^{\rm thin}\gg\delta t_{\rm rad}^{\rm thick}$, so their sum is determined by $\delta t_{\rm rad}^{\rm thin}$, as expected; see figure 5 for small values of $\ln p$, where the blue line for $\delta t_{\rm rad}^{\rm thin}$ is the highest. In the deeper layers below the photosphere, the situation is the other way around and the time step is expected to be governed by $\delta t_{\rm rad}^{\rm thick}$, again as expected; see figure 5 for large values of $\ln p$ (i.e., below the photosphere), where the red line for $\delta t_{\rm rad}^{\rm thick}$ is the highest. The radiative time step is thus determined by the minimum of $\delta t_{\rm rad}$, corresponding to the black dashed line in figure 5. We find that $\delta t_{\rm rad}$ is the shortest either in the photosphere [panels (a) and (b) of figure 5], indicated by a red dot where $\tau=1$, or just below it [see panels (c) and (d) of figure 5, i.e., for models with shallow or no outer convection zone].

Note that our restriction of static diffusion may no longer be justified if $\beta\tau\gg 1$ in the deeper layer of the stars, which are known to be in the dynamic diffusion regime; see Krumholz et al. (2007). To estimate the depth where this occurs in the stars in figure 5, let us estimate $u$ based on the assumption that the convective flux can be modelled by the associated mixing length expression, $\sigma_{\rm SB}T_{\rm eff}^{4}=\rho u^{3}$. The depth below which dynamic diffusion occurs is obtained when $\beta\tau\approx(\sigma_{\rm SB}T_{\rm eff}^{4}/\rho)^{1/3}\tau/c$ exceeds unity. In that regime, our radiation transport equations are no longer applicable, which is denoted by the grey regions in figure 5. We see that $\beta\tau\ll 1$ is always obeyed in the surface layers of all the stars, where our radiation transport equations are valid. Note that our estimate of $\beta\tau$ applies only to regions where convection is possible, i.e., in the outer layers of cool stars or in the convective cores of hot stars. Thus, the actual depth where $\beta\tau\approx 1$ may differ.

The maximum permissible time step $\delta t$ for this problem is then given by the minimum of $\delta t_{\rm rad}$ and the acoustic time step $\delta t_{\rm s}$ (the green solid line). As can be seen from figure 5, for $T_{\rm eff}=5000\,{\rm K}$, it is always $\delta t_{\rm s}$ that constrains the actual time step. Although for $T_{\rm eff}=7000\,{\rm K}$ the minimum of $\delta t_{\rm rad}$ occurs below the photosphere, at that location $\delta t_{\rm s}$ is still the shorter time step. However, for $T_{\rm eff}=10,000\,{\rm K}$ and higher, $\delta t_{\rm s}$ is no longer the constraining time step, instead it is $\min(\delta t_{\rm rad})$. To better understand the dependence of various time steps on $T_{\rm eff}$ and $g$, we turn to table LABEL:THR and figure 6. In table LABEL:THR, we give the radiative and acoustic time steps in the photosphere, $\delta t_{\rm rad}^{\rm phot}$ and $\delta t_{\rm s}^{\rm phot}$, respectively, and the radiative and acoustic time steps at the position where the radiative time step has a minimum, $\delta t_{\rm rad}^{\rm deep}$ and $\delta t_{\rm s}^{\rm deep}$, respectively; for the five models (i)–(v) discussed above in figures 4 and 5, and for a few other ones. In figure 6, we plot these newly defined time steps as a function of $T_{\rm eff}$ for different values of $g$.

We see from table LABEL:THR and figure 6 clearly that as $T_{\rm eff}$ increases for a given $g$, the radiative time steps (thick lines) decrease whereas the acoustic time steps (thin lines) increase (both in the photosphere and deeper layers). We find that for stars with $T_{\rm eff}<10,000\,{\rm K}$, $\delta t_{\rm s}$ gives the time step constraint, whereas for stars with $T_{\rm eff}\gtrsim 10,000\,{\rm K}$, $\delta t_{\rm rad}$ is the more constraining one. Note that this segregation of the cold and hot stellar branches at $T_{\rm eff}\approx 10,000\,{\rm K}$ seen in the $\delta t_{\rm rad}$ curves of figure 6 is a consequence of our opacity prescription given by (28). In this equation, we use $T_{0}=13,000\,{\rm K}$, such that the stars hotter than this are described predominantly by Kramers opacity and stars cooler than this by H- opacity. We see from figure 6 that for a given $g$, the $\delta t_{\rm s}$ and $\delta t_{\rm rad}$ curves intersect at a particular $T_{\rm eff}$, which is typically $<10,000\,{\rm K}$; for e.g. the thin and thick red dashed lines in panel (a) intersect at $\sim 9000\,{\rm K}$. These stars are in fact the most economical to simulate numerically, as $\delta t_{\rm s}\approx\delta t_{\rm rad}$ and one does not need to worry about conflicting time steps in the problem. For hotter stars, however, we see that the $\delta t_{\rm s}$ and $\delta t_{\rm rad}$ diverge more and more away from each other, thus leading to a problem. For a given $T_{\rm eff}$, we see that $\delta t_{\rm s}$ always decreases with increasing $g$. The $\delta t_{\rm rad}$ curves for a given $T_{\rm eff}<10,000\,{\rm K}$, on the other hand, seem to be nearly independent of $g$. However, as $T_{\rm eff}$ increases beyond $10,000\,{\rm K}$, we find that the more massive stars have a shorter radiative time step. Interestingly, as $T_{\rm eff}\rightarrow 10^{5}\,{\rm K}$, $\delta t_{\rm rad}$ again starts to become independent of $g$. Note that the time steps in the deeper layers of the star, as shown in figure 6(b) are much shorter than their photospheric counterparts, especially for hot massive stars. This also poses a numerical challenge if one wishes to simulate such a star all the way from the deeper layers up to the photosphere. Figure 6 can be thought to represent a Hertzsprung–Russell (HR) diagram for stars having different $g$ and $T_{\rm eff}$.

Finally, we confirm that the stratification in the present hydrostatic models agrees with the steady state solution obtained in section 4.1 using the Pencil Code. We do so by comparing the locations of the photosphere in figures 3 and 5(d), both of which correspond to the model having $T_{\rm eff}\sim 70,000\,{\rm K}$. We note from figure 3 that the computational domain ranges from $z=-30\,{\rm Mm}$ to $+30\,{\rm Mm}$, with the negative values representing the deeper layers of the star. The corresponding $\log_{10}p$ varies from $2.75$ to $5.47$ across the domain, where $p$ is in cgs units. On comparing with figure 5(d), we see that this is well within the regime of static diffusion, where our numerical calculations are valid. In figure 3, the photosphere ($\tau=1$) is at $z=-5.8\,{\rm Mm}$ or $\log_{10}p=4.6$, which indeed is in excellent agreement with the location of the photosphere in figure 5(d). Also, the value of $\delta t_{\rm rad}$ at the photosphere agrees between figures 3 and 5(d), being $\approx 10^{-2}\,{\rm s}$ in both (keeping in mind that the time step in figure 3 is given in milliseconds). We make a note here regarding the location of the minimum of the radiative time step in theoretical models versus numerical computation. We see from figure 5(d) that $\min(\delta t_{\rm rad})$ occurs in the deeper layers at $\log_{10}p\approx 7.2$, which is in fact outside the computational domain of interest in figure 3. It is only for the cool stars, that the photosphere coincides with the location of $\min(\delta t_{\rm rad})$; see e.g. figure 5(a).

4.3 Simulation results for disc models

In this section, we solve for the structure of an accretion disc around a white dwarf of mass $M_{\rm WD}=1.1M_{\odot}$. We implement the shearing box method in the Pencil code, such that the disc is located at a distance of $r_{\rm disc}=10^{10}\,{\rm cm}$ from the central star, with $\varOmega=0.012~{}{\rm rad}\,{\rm s}^{-1}$ at this location. Thus, one rotation period at $r_{\rm disc}$ is $T_{\rm rot}=2\pi/\varOmega=0.52\,{\rm ks}$. In our earlier exploratory models, we solved the disc models for both the upper and lower disc plane, but in the subsequent models presented below, we have restricted ourselves to solving the equations in just the upper disc plane by assuming a symmetry condition at $z=0$. These solutions are identical to our earlier ones and computationally more economic.

We consider one-dimensional solutions of equations (13)–(16) for a fixed value of the vertically integrated (surface) density $\Sigma=\int_{-\infty}^{\infty}\rho\,{\rm d}{}z$, with $\nu=5\times 10^{10}\,{\rm cm}^{2}\,{\rm s}^{-1}$ and $C_{\rm shock}=3$. In our horizontally periodic domains with horizontal extent $L_{x}\times L_{y}$, the total mass is $\Sigma L_{x}L_{y}$, which is conserved and therefore given by the initial conditions. Again, see appendix A and appendix C for details about the boundary conditions on various variables. We use $\Sigma=7\times 10^{-6}\,{\rm g}\,{\rm cm}^{-3}\,{\rm Mm}$, and two values for ${\cal H}_{0}={\cal H}(z=0)$, namely, ${\cal H}_{0}=2\times 10^{-6}\,{\rm g}\,{\rm cm}^{-3}\,{\rm km}^{3}\,{\rm s}^{-3}\,{\rm Mm}^{-1}$ (cold disc) and ${\cal H}_{0}=5\times 10^{-5}\,{\rm g}\,{\rm cm}^{-3}\,{\rm km}^{3}\,{\rm s}^{-3}\,{\rm Mm}^{-1}$ (hot disc) and $z_{\rm heat}=1\,{\rm Mm}$ for both cases. For the initial isothermal stratification used here (as discussed in section 3.5), the pressure and density scale heights are identical, so $H_{\!p}=c_{\rm s0}/\varOmega=H_{\!\rho}$, but in general they are somewhat different from each other; see appendix D. For the cold disc, we set $H_{\!p}=0.42\,{\rm Mm}$ ($c_{\rm s0}=5\,{\rm km}\,{\rm s}^{-1}$) and for the hot disc, $H_{\!p}=2.5\,{\rm Mm}$ ($c_{\rm s0}=30\,{\rm km}\,{\rm s}^{-1}$). In cgs units, our values correspond to $\Sigma=7\times 10^{2}\,{\rm g}\,{\rm cm}^{-2}$, ${\cal H}_{0}=20\,{\rm g}\,{\rm cm}^{-1}\,{\rm s}^{-3}$, and ${\cal H}_{0}=500\,{\rm g}\,{\rm cm}^{-1}\,{\rm s}^{-3}$, which are appropriate values for discs in cataclysmic variables. The corresponding accretion rates can be obtained from (29), and are given by $\dot{M}=5.8\times 10^{13}\,{\rm g}\,{\rm s}^{-1}$ (cold disc) and $\dot{M}=1.4\times 10^{15}\,{\rm g}\,{\rm s}^{-1}$ (hot disc). These were chosen based on preliminary calculations of semi-analytically constructed models relevant to the regime where multi-valued solutions of $\dot{M}$ (or equivalently ${\cal H}_{0}$ in our formalism) are possible for a given $\Sigma$. These solutions are governed by the hydrogen ionisation instability that lead to the so called disc instability model of cataclysmic variables; see Lasota (2001) for a review. Note that for Pencil Code simulations, it is convenient to measure lengths in units of $\,{\rm Mm}$, speed in $\,{\rm km}\,{\rm s}^{-1}$, density in $\,{\rm g}\,{\rm cm}^{-3}$, time in $\,{\rm ks}$, and temperature in Kelvin, which explains our choice of units adopted in the results presented here. Other choices would have been possible, too, although working with $\,{\rm Mm}$ and $\,{\rm km}\,{\rm s}^{-1}$ is useful because those are the units displayed in many of the diagrams.

We first discuss figure 7, where we present the solutions for a cold accretion disc. In this case, $T$ varies from a photospheric value of about $2600\,{\rm K}$ to about $3600\,{\rm K}$ in the midplane; see panel (a). Here we choose a vertical domain size of $L_{z}=1\,{\rm Mm}$. The photosphere (i.e., $\tau=1$ surface), as denoted by the red dots in the figure, occurs at a depth of $0.65\,{\rm Mm}$. In panels (c) and (d), we plot the various time steps for two vertical grid resolutions $N_{z}=144$ and $576$, which were run for a total time of $3\,{\rm ks}$ ($5.7~{}T_{\rm rot}$) and $5\,{\rm ks}$ ($9.7~{}T_{\rm rot}$), respectively. According to our reasoning in section 2, the radiative time step should be limited by the sum of $\delta t_{\rm rad}^{\rm thin}$ and $\delta t_{\rm rad}^{\rm thick}$. However, for the cold disc solution, the former is much larger. Here, $\delta t_{\rm rad}^{\rm thin}$ reaches values of between $10^{-2}\,{\rm ks}$ in the midplane and increases to $10^{-1}\,{\rm ks}$ in the outer layers (for both $N_{z}=144$ and $576$ as it is independent of $\delta z$), and, hence it determines the radiative time step $\delta t_{\rm rad}$. The maximum permissible empirically determined time step $\delta t$, on the other hand, is only around $10^{-3}\,{\rm ks}$ for $N_{z}=144$, and $1.8\times 10^{-4}\,{\rm ks}$ for $N_{z}=576$. It turns out that this value is entirely explained by the standard CFL condition, where the time step is limited by the acoustic time step $\delta t_{\rm s}$; compare black dot-dashed and green solid lines in figures 7(c) and (d). From panel (b) we see that $c_{\gamma}$ is well below $c_{\rm s}$ throughout the computational domain. It varies only between $1.6\,{\rm km}\,{\rm s}^{-1}$ in the photosphere to about $2.6\,{\rm km}\,{\rm s}^{-1}$ in the midplane, and reaches $6\,{\rm km}\,{\rm s}^{-1}$ in the outermost parts well above the photosphere. The relevance of the ratio $c_{\rm s}/c_{\gamma}$ will become more clear in section 5.

In figure 8, we present the hot disc solutions, where $T$ varies from a photospheric value of $5000\,{\rm K}$ to about $32,000\,{\rm K}$ in the midplane. We choose $L_{z}=4\,{\rm Mm}$ for this case, which extends from the midplane to a point somewhat above the photosphere. The photosphere in this case occurs at a depth of $2.8\,{\rm Mm}$. We again show the results for $N_{z}=144$ and $576$, which were run for a total duration of $4\,{\rm ks}$ ($7.7~{}T_{\rm rot}$) and $0.018\,{\rm ks}$ ($0.03~{}T_{\rm rot}$) respectively. We find from panels (c) and (d) that the hot disc behaves somewhat similarly to the stellar surface layers discussed in section 4.2, where $\delta t_{\rm rad}$ in the layers at some depth beneath the photosphere is governed by $\delta t_{\rm rad}^{\rm thick}$ (solid red line), while in the outer layers by $\delta t_{\rm rad}^{\rm thin}$ (solid blue line). The minimum of the total radiative time step, however, occurs at the midplane, such that $\min(\delta t_{\rm rad})\sim 10^{-4}\,{\rm ks}$ for $N_{z}=144$, and $\sim 7\times 10^{-6}\,{\rm ks}$ for $N_{z}=576$. This is very much consistent with the maximum permissible empirical time step $\delta t$, which is $\sim 6\times 10^{-5}\,{\rm ks}$ for $N_{z}=144$ and $\sim 7\times 10^{-6}\,{\rm ks}$ for $N_{z}=576$. Note that $\delta t_{\rm rad}^{\rm thick}$ is much more stringent in this case and limits the time step; compare the red solid and black dot-dashed lines in figure 8(d). Also, $\delta t_{\rm rad}^{\rm thick}\propto\delta z^{2}$, as can be seen from panels (c) and (d) of both figures 7 and 8. The standard CFL condition is not relevant to explain $\delta t$ in this case as $\delta t_{\rm s}>\min(\delta t_{\rm rad})$. Furthermore, $c_{\gamma}>c_{s}$ throughout the domain, reaching about $7000\,{\rm km}\,{\rm s}^{-1}$ in the midplane and even $10^{5}\,{\rm km}\,{\rm s}^{-1}$ in the outer layers; see figure 8(b). We will return to the effect of this in section 5.

Finally, we discuss table LABEL:Tsum, where we summarise the empirically determined maximally permissible time step $\delta t$ for the cold and hot disc models. We also compare with the numerically determined values of $\chi_{0}\delta t/\delta z^{2}$ and $c_{\rm s0}\delta t/\delta z$, which will help us constrain the coefficients $C_{\rm rad}^{\rm thick}$ and $C_{\rm CFL}$, respectively (the subscript 0 indicates the values of the respective quantities at the disc midplane). From the cold disc solutions, where the time step is constrained by $\delta t_{\rm s}$, we find that $c_{\rm s0}\delta t/\delta z=1.3$ for $N_{z}=144$, and $c_{\rm s0}\delta t/\delta z=0.95$ for $N_{z}=576$. These values are very close to the value of the standard Courant factor used in the Pencil Code. Hence, we can conclude that $C_{\rm CFL}=0.95$, which is also consistent with the value of 0.9, which was adopted while plotting figure 7. Note that ideally these values should be independent of resolution, but in practice they do seem to depend on it. We hence adopt the smaller of the two values as the more restrictive constraint. From the hot disc solutions, where the time step is constrained by $\delta t_{\rm rad}^{\rm thick}$, we hope to constrain the coefficient $C_{\rm rad}^{\rm thick}$. First we find that for this case, $c_{\rm s0}\delta t/\delta z\ll 1$, which is contrary to our understanding of the CFL coefficient and, hence, we discard these values. Comparing the values of $\chi_{0}\delta t/\delta z^{2}$, we conclude that $C_{\rm rad}^{\rm thick}=0.19$, since this is the smallest of the two resolutions. This value is also consistent with our choice of 0.2, which was used while plotting figure 8. Note that neither the hot disc nor the cold disc models are suitable for determining the coefficient $C_{\rm rad}^{\rm thin}$, as $\delta t_{\rm rad}^{\rm thin}$ never constrains the time step in these cases. However, the value of $C_{\rm rad}^{\rm thin}=4$ that we obtained from our earlier one-dimensional experiment and used for our plots, is indeed consistent with the simulation results reported in section 4.1; see also figure 3.

4.4 Three-dimensional simulations of discs with vertical heating profile

In section 3.5, we introduced the heating profile ${\cal H}(z)$ and stated that the discontinuities at $z=\pm z_{\rm heat}$ do not cause any artefacts in the thermodynamic variables such as temperature and specific entropy. This is demonstrated in figure 9, where we show the results of a three-dimensional simulation (of both the gas and the radiation field) for the same white dwarf-accretion disc system as in section 4.3, but with an intermediate heating source ${\cal H}_{0}=10^{-5}\,{\rm g}\,{\rm cm}^{-3}\,{\rm km}^{3}\,{\rm s}^{-3}\,{\rm Mm}^{-1}$ (or equivalently, $\dot{M}=2.9\times 10^{14}\,{\rm g}\,{\rm s}^{-1}$) in a cubic domain $[L_{x},L_{y},L_{z}]$ of size $(4\,{\rm Mm})^{3}$, covering both the lower and upper disc planes, using $576^{3}$ mesh points, $z_{\rm heat}=1\,{\rm Mm}$, $\nu=5\times 10^{9}\,{\rm cm}^{2}\,{\rm s}^{-1}$ and $C_{\rm shock}=0.5$. Again, see appendices A and C for details about the boundary conditions on various variables. The 3D simulation was run for a total time of $27.5\,{\rm ks}$ or $52.8~{}T_{\rm rot}$.

The temperature reaches a maximum of about $8000\,{\rm K}$ at $z=0$ and has an approximately flat profile away from the lower and upper photospheres (i.e., at $\pm 1.3\,{\rm Mm}$). It turns out that the $z$ dependence of temperature profile is not perfectly flat. This can probably be ascribed to the long thermal adjustment time in this system, which we can estimate as follows. The Kelvin–Helmholtz timescale is given by $\tau_{\rm KH}=E_{\rm th}/L$, where $E_{\rm th}=\int\rho c_{\!v}T\,{\rm d}{}V$ is the internal energy, $L=2F_{\infty}L_{x}L_{y}=2{\cal H}V$ is the luminosity for the losses on both photospheres, $V=L_{x}L_{y}L_{z}$ is the volume of the domain, and $F_{\infty}$ is the value of the total flux in the photosphere. Using $E_{\rm th}=\Sigma L_{x}L_{y}\overline{c_{\rm s}^{2}}/[\gamma(\gamma-1)]$ with $\overline{c_{\rm s}^{2}}\equiv\langle\rho c_{\rm s}^{2}\rangle/\langle\rho\rangle\approx(13\,{\rm km}/\,{\rm s})^{2}$, and $[\gamma(\gamma-1)]^{-1}=0.9$ for $\gamma=5/3$, we have $\tau_{\rm KH}=0.45\,\overline{c_{\rm s}^{2}}\Sigma/{\cal H}_{0}L_{z}=50\,{\rm ks}$, which is about three times longer than the duration of our simulation.

The specific entropy, defined here as $s/c_{\!p}=\ln(T/p^{\nabla_{\rm ad}})$, has a negative slope, ${\rm d}{}s/{\rm d}{}\ln z<0$, corresponding to a Schwarzschild-unstable stratification. This should lead to instability and hence to turbulent convection. This is indeed the case as seen in figure 9(c), where we plot the mean energy fluxes averaged over a time span when the system is in an approximately steady state between $t=15\,{\rm ks}$ and $20\,{\rm ks}$, which is still subject to slow thermodynamic adjustments on the longer timescale $\tau_{\rm KH}=50\,{\rm ks}$. In particular, we plot $F_{\rm rad}$ together with $F_{\rm conv}$ and $F_{\rm tot}=F_{\rm rad}+F_{\rm conv}$, where $F_{\rm conv}=F_{\rm enth}+F_{\rm kin}$, with

[TABLE]

being the enthalpy and kinetic energy fluxes, and overbars denote from now on horizontal $xy$ averages.

Note that in figure 9, $F_{\rm tot}$ varies approximately linearly in $|z|<z_{\rm heat}$ and then reaches a plateau with $F_{\rm tot}=\pm F_{\rm tot}^{\infty}$ on both ends. However, $F_{\rm tot}^{\infty}/{\cal H}_{0}z_{\rm heat}\approx 1.2$ exceeds the expected values of $\pm 1$, which is, again, indicative of the simulation still not being in thermal equilibrium. The radiative flux is small near the midplane and energy is mostly carried by convection. The kinetic energy is not plotted, but it is about 10% of the convective flux and directed opposite to it, i.e., inward. This is a well known consequence of the up-down asymmetry of compressible convection (Hurlburt et al., 1984).

We see that in the inner parts of the disc, i.e., for $|z|<0.5\,{\rm Mm}$, most of the energy is carried by convection. This is curious in view of the fact that the stratification in those parts is close to adiabatic. Therefore, the standard mixing length prescription of the convective flux being carried by a gradient term, $F_{\rm enth}\approx F_{\rm G}$ with (e.g. Rüdiger, 1989)

[TABLE]

where $\chi_{\rm turb}\approx u_{\rm rms}H_{\!p}/3$, cannot hold. As argued in Brandenburg (2016), the reason for this is that there is another important term, the Deardorff term, resulting from entropy fluctuations (Deardorff, 1966, 1972),

[TABLE]

where we have ignored the possibility of factors of the order of unity. Likewise, in (33), there could also be such factors, so we cannot expect perfect agreement between $F_{\rm conv}$ and the contributions $F_{\rm G}$ and $F_{\rm D}$ shown in figure 9(c). We do see, however, that $F_{\rm D}$ increases approximately linearly in the bulk of the disc ($|z|<0.5\,{\rm Mm}$), while $F_{\rm G}$ does not even have the correct sign in order to explain $F_{\rm conv}$. This is strong evidence that convection in such discs must be described by the Deardorff term. Similar results have previously only been found from stellar convection simulations (Käpylä et al., 2017, 2019).

Returning to the topic of this paper, we now investigate the various times step constraints for this model. In figure 10, similar to our earlier plots, we show the various times steps, but now we also plot the time step constraint from viscosity, particularly the one from shock viscosity. Note that the empirical time step $\delta t$ in this case is not fixed at the beginning of the simulation (as done in figures 3, 7 and 8), but is allowed to evolve. It turns out that the minimum of radiative and acoustic time steps are approximately equal ($\approx 5\times 10^{-4}\,{\rm ks}$) and occur in the disc midplane ($z=0$). However, the most severe time step constraint comes, in this case, from the outer layers above the photosphere, where the shock viscosity is large and hence the corresponding time step is small (green dashed line). This is in accord with the empirical time step $\delta t\approx 8\times 10^{-5}\,{\rm ks}$ attained at the end of the the simulation.

To appreciate the reason for such shocks to occur, we now look at the $xz$ cross sections of $u_{x}$, $u_{y}$, and $u_{z}$ shown in figure 11. We see that the flow possesses two major pairs of up- and down-draughts. These cross sections were taken at $y=0$, but different cross sections for different values of $y$ look qualitatively similar, indicating that the large-scale structure is independent of $y$. Away from the both photospheres, however, significant velocities are still possible, which can lead to the formation of shocks in those parts. These are best seen in the image of $u_{x}$; see first panel of figure 11. The arrow points to the strongest shock near $z=-2\,{\rm Mm}$, where the local Mach number, i.e., the ratio of velocity to sound speed, reaches values of around 1.3.

5 Approaches to mitigating the radiative time step problem

There appears to be a general problem with hot stars and accretion discs where the radiative time step is much shorter than the acoustic one. This is either because

[TABLE]

as in section 4.1, or because

[TABLE]

as in section 4.3. In both cases, the problem is caused by the smallness of the factor $\mbox{\rm Bo}=16c_{\rm s}/c_{\gamma}$ compared with the value of $\delta z/4\ell$ in problem A or the value of $\ell/0.6\delta z$ in problem B. In problem A, $c_{\rm s}/c_{\gamma}$ becomes as small as $10^{-6}$, as seen by comparing the green solid and dotted lines at the outer edge of the computational domain in figure 3(b). While in problem B, $c_{\rm s}/c_{\gamma}$ drops to $10^{-4}$ at the outer edge of the disc, as seen in figure 8(b). For problem A, the photospheric value of $\delta z/4\ell$ is about 0.006, as calculated from the model discussed in figure 3. While for problem B, the photospheric value of $\ell/0.6\delta z$ is about $16$, as calculated from the model discussed in figure 8(d). Note that these are both quite large compared to the respective $c_{\rm s}/c_{\gamma}$ values.

To discuss ways of mitigating the time step problem, we must distinguish the two cases A and B. For problem B, the difficulty appears through the optically thick constraint in the part where convection would develop in the three-dimensional simulations, as discussed in section 4.4; see also figure 11. Thus, radiation should not be physically important and it would be unreasonable to spend much computational resources on this. One may therefore adopt the approach of artificially increasing the opacity in that part, so that $\chi=c_{\gamma}\ell/3=c_{\gamma}/3\kappa\rho$ becomes smaller and, as a consequence, the optically thick time step becomes longer. This would increase the fractional convective flux, but not the total flux.

For problem A, on the other hand, it is important to maintain a large radiation pressure in order to study photoconvection and, hence, $\kappa F_{\rm rad}/cg$ should be close to unity. Here, however, increasing $\kappa$ would be counterproductive, because the relevant time step is $4\,\ell/c_{\gamma}$, which would become even smaller as $l\propto 1/\kappa$. Thus, we have the problem that, on the one hand, $\kappa F_{\rm rad}/cg$ should be close to unity and certainly not be too small, and on the other hand, $4\,\ell\,c_{\rm s}/(0.9\,\delta z\,c_{\gamma})$ should also be close to unity so that the acoustic and optically thin radiative time steps are close together; see equation (35). Combining these two constraints, we want to ensure that the

[TABLE]

is as close to unity as possible because we want $\kappa F_{\rm rad}/cg\gg 1$. Using $F_{\rm rad}\approx\sigma_{\rm SB}T^{4}$, $\delta z\approx 0.05H_{\!p}$, and $\kappa\rho\ell=1$, we find for the

[TABLE]

Thus, it is clear that the best chance of studying photoconvection in hot massive stars is in the relativistic regime, because then the right hand side of equation (38) is close to unity. In that case, however, one can no longer neglect the time dependence of the radiation field. Alternatively, of course, one may just take a solar-type model, where the radiative and acoustic time steps are already close together, but then decrease $c$ artificially, in order to boost the radiative pressure and reduce computational costs. This approach is similar to the reduced speed of light approximation (RSLA) originally proposed by Gnedin and Abel (2001). It has been employed, for e.g., by Skinner and Ostriker (2013), who adopt a particular closure relation and solve the radiation momentum equation semi-explicitly, in a module for the Athena code (Stone et al., 2008). They found that the RSLA method is best applicable for systems with moderate optical depth ($\lesssim 10$), such as several star formation regimes in galactic discs.

At the end of section 2, we mentioned the possibility that also a strong radiation pressure could in principle restrict the time step. Such a situation may arise, for example, in the deep interior of a super-massive star; see figure 5. A relevant application of the RSLA is therefore the study of such a constraint. We therefore now reconsider the simulation of the hot star of section 4.1, where the radiation pressure was retained in the momentum equation (14). We artificially lower the value of $c$ in the term $g_{\rm rad}=\kappa F_{\rm rad}/c$ of this equation to make it more pronounced. We find that it is only when we lower it from $c=3\times 10^{5}\,{\rm km}\,{\rm s}^{-1}$ to about $0.02\,{\rm km}\,{\rm s}^{-1}$ that this term begins to restrict the time step.

In table 4, we list such models, where we lower the value of $c$ and increase $g$ such that $g\approx g_{\rm rad}$ so that the model remains in Eddington equilibrium, i.e., $g_{\rm eff}\to 0$ in equation (21). Note that all the other parameters of the simulation remain the same as mentioned in section 4.1. The time needed to accelerate a parcel of gas over a distance $\delta z$ is $\sqrt{2\delta z/g_{\rm rad}}$. We can therefore formulate a new time step constraint by generalising (11) to

[TABLE]

with a new coefficient $C_{\rm rad}^{\rm pres}$ and the ellipsis indicating the other terms of equation (11). This additional constraint is qualitatively different from any of our earlier constraints in that it is proportional to $\delta z^{1/2}$. The results presented in table 4 suggest that $C_{\rm rad}^{\rm pres}\approx 0.4$. This new finding illustrates the potential usefulness of applying the RSLA.

By changing or modifying the physical problem, one is clearly never able to solve the actual problem, but, under certain circumstances, one may be able to address relevant aspects of the problem. An example of this type of approach is the work of Käpylä et al. (2013, 2019), who study astrophysically relevant aspects of rotating convection by firstly increasing the radiative flux to bring the Kelvin-Helmholtz timescale closer to other timescales in the problem, and secondly increasing the angular velocity such that the rotational influence on the flow is in the relevant regime. How useful such approaches are cannot easily be answered and may ultimately need to await a full solution to the problem. Another example is the study of high Rayleigh and Reynolds number flows, where we will never be able to solve the astrophysically relevant regime using brute force approach.

6 Conclusions

Our work has demonstrated quantitatively how the maximally permissible time step varies with $\delta z^{2}/\chi D$ in the optically thick regime and with $\ell/c_{\gamma}$, independent of $\delta z$, in the optically thin regime. In particular, we have shown that the radiative time step is governed by the sum of the optically thick and thin constraints, and not by the smaller of the two, as one might have naively expected by analogy with other time step constraints entering the problem. This agrees with the work of Freytag et al. (2012), except that our studies identify the existence of two independent coefficients in front of the terms for the optically thin and thick regimes. It should be noted, however, that the large difference between the coefficients $C_{\rm rad}^{\rm thick}$ and $C_{\rm rad}^{\rm thin}$ could have been alleviated by working with the Nyquist wavenumber instead of the inverse mesh spacing. This amounts to including a $\pi^{2}$ factor in the definition of $C_{\rm rad}^{\rm thick}$; see our discussion in section 2. However, to facilitate comparison with earlier work, we have refrained from redefining this coefficient.

The fact that we find $C_{\rm rad}^{\rm thick}\approx 0.2$ and $C_{\rm rad}^{\rm thin}\approx 4$ over a broad range of different physical circumstances supports our expectation that these values are universal. While this is to be expected, we should nevertheless keep an open mind and continue verifying this in new and as yet unexplored regimes whenever possible.

We have also demonstrated that a radiative time step constraint of the form $\delta z/c_{\gamma}$, analogous to the Courant condition as given by (10), is not, in general, justified; see especially our figures 3 and 8. We thank Matthias Steffen (private communication) for pointing out to us that such a relation could be justified in principle by noting that (4) predicts maximum cooling when $\ell=\sqrt{3}/k_{\rm Ny}\approx 0.55\,\delta z$, so that the cooling time becomes $2\ell/c_{\gamma}\approx 1.1\delta z/c_{\gamma}$. The simulation presented in figure 3 has demonstrated, however, that this would result in a time step constraint that would be about 1000 times shorter (see minimum of green dotted line in figure 3) than the actual maximum permissible time step (see black dot-dashed line) required for stability. Thus, $\delta z/c_{\gamma}$ does not appear to be a useful indicator of the maximum permissible time step.

Incidently, figure 5(a) shows that, in the proximity of the photosphere, the local value of $\delta z/c_{\gamma}$ is indeed close to the actual maximum permissible time step of about $10^{4}\,{\rm s}$. However, already in figure 5(b) this is no longer the case, because the black dashed line attains its minimum well below the photosphere, where $\log_{10}p\approx 8$. Furthermore, in figure 5(a), the expected time step constraint is apparently not given by the local minimum of $\delta z/c_{\gamma}$, but rather by the local maximum of $\delta z/c_{\gamma}$. This is here because $c_{\gamma}$ has a minimum close to the photosphere. Deeper down it increases because $T$ increases, and higher up it also increases because $\rho$ drops rapidly.

By covering stellar surface models in the HR diagram, we get a comprehensive understanding of which of the two constraints decide about the limiting time step for a large range of different circumstances. This showed that the most severe constraint on the radiative time step occurs for larger values of $g$ and $T_{\rm eff}$. This corresponds to the lower left corner of the HR diagram, recalling that in theoretical HR diagrams, increasing luminosity is replaced by decreasing surface gravity; see, e.g., Trampedach et al. (2013). We have also seen that, for the cool stars with $T_{\rm eff}\lesssim 5000\,{\rm K}$ the location of the minimum radiative time step coincides with the location of the photosphere. This is not true for hotter stars, where the minimum radiative time step occurs in the deeper layers. For accretion discs (both hot and cold), on the other hand, the shortest radiative time step tends to occur in the midplane.

The examples presented in this paper highlight some of the difficulties in dealing with global simulations by covering regimes where different time step constraints prevail. It is clear that the optimal approach would be one where different regions in space would not only have different spatial resolutions, but also different time steps. This would save resources that can at the same time be used to speed up the calculation in regions that require shorter time steps. A code satisfying such requirements is the DISPATCH code (Nordlund et al., 2018). Nevertheless, the time step constraints discussed in the present paper should apply to such codes just as well.

For the three dimensional accretion disc model presented here, the temperatures are moderate, and the acoustic and radiative time step constraints are about equally short. This happens in the midplane. In the outer parts above the two photospheres, shocks become important. The time step constraint resulting from the shock viscosity is here about equally severe as the acoustic and radiative ones in the disc midplane. In such a situation, one might not gain much by using an implicit scheme for radiation, or by adopting a code that treats different regions in space with different times steps. However, it is important to monitor the various times step constraints carefully and try to stay close to physical regimes in which the different constraints are not vastly different from each other. In this way, the simulation can utilise existing resources in an optimal way.

Acknowledgements

This paper is dedicated to Ed Spiegel. If it was not for the radiative time step problem, the Brandenburg and Spiegel (2006) paper would have been published by now! We thank Matthias Steffen for sharing with us his experience and knowledge regarding the time step constraint in stellar convection simulations. We are also grateful to the three referees for their thoughtful comments. AB also acknowledges Fazeleh (Sepideh) Khajenabi for her work on the radiative time step problem while visiting Nordita in the spring of 2010. UD thanks Akshay Bhatnagar for useful discussions. This research was supported in part by the National Science Foundation under the Astronomy and Astrophysics Grants Program (grant 1615100), and the University of Colorado through its support of the George Ellery Hale visiting faculty appointment. Simulations presented in this work have been performed with computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm. The source code used for the simulations of this study, the Pencil Code, is freely available on https://github.com/pencil-code/. The DOI of the code is http://doi.org/10.5281/zenodo.2315093. The setups of runs and corresponding data are freely available on https://www.nordita.org/~brandenb/projects/tstep/.

Appendix A Radiative boundary conditions

In section 3.2, we stated that the assumption of zero incoming intensity is not accurate when we want to reproduce the analytic solution for an infinitely extended layer, where the gas beyond the simulated boundary does contribute to producing incoming radiation. To take this into account, we assume that $F_{\rm rad}$ is known and that the system is in radiative equilibrium. In that case, and with just two rays, I_{\pm}=I(\bm{x},t,\pm\hat{\mbox{\boldmathz}}{}), we have

[TABLE]

so $I_{\pm}=(S\pm F_{\rm rad}/4\pi)/2$. In the Pencil Code, this boundary condition is invoked by stating the symbolic name bc_rad=‘p’,‘p’,‘S+F:S-F’. The three entries correspond to the three directions $(x,y,z)$ respectively, ‘p’ represents periodic boundary condition, S+F applies to the lower $z$ boundary and S-F applies to the upper $z$ boundary (and are thus separated by a colon). Note that these are the boundary conditions used for the hot stellar model in figure 3.

The radiative boundary conditions used for the 1D accretion disc models in figures 7 and 8 are bc_rad=‘p’,‘p’,‘S:0’, where 0 represents zero value in ghost zones and free value on the upper $z$ boundary. The conditions used for the three dimensional accretion disc model in section 4.4 are bc_rad=‘p’,‘p’,‘0’.

Appendix B Different formulations of the energy equation

The energy equation can be formulated in a number of equivalent forms. Equation (15) illustrates that heating and cooling only affect the specific entropy, while equation (3) illustrates that temperature change is governed by the specific heat at constant pressure when ${\rm D}{}p/{\rm D}{}t$ can be neglected, which is when sound waves equilibrate pressure fluctuations. Another formulation of the right-hand side of equation (15) is

[TABLE]

which shows that temperature changes depend on the specific heat at constant volume when ${\rm D}{}\ln\rho/{\rm D}{}t$ can be neglected. The formulation in (3) has been favoured by Spiegel and Veronis (1960) to show that in the Boussinesq approximation, where \mbox{\boldmath\nabla}{}\bm{\cdot}\bm{u}=0, the relation {\rm D}{}\ln\rho/{\rm D}{}t=-\mbox{\boldmath\nabla}{}\bm{\cdot}\bm{u} is no longer invoked and ${\rm D}{}\ln\rho/{\rm D}{}t$ cannot be neglected. It is this fact that also motivates the presence of the $\gamma=c_{\!p}/c_{\!v}$ factor in the definition of the Prandtl number $\Pr=\nu/\gamma\chi$.

Appendix C Boundary conditions for hydrodynamic variables

The purpose of this appendix is to discuss numerical details regarding the models presented in section 4.1, section 4.3 and section 4.4.

The independent variables solved for the hot stellar surface model in figure 3 are $u_{x}$, $u_{y}$, $u_{z}$, $\ln\rho$, $\ln T$ and $\nu_{\rm shock}$. The boundary conditions supplied for these 6 variables along the $z$-direction is given in Pencil Code notation as, bcz=‘s’,‘s’,‘a’,‘e2’,‘e2’,‘s’; where ‘s’ implies symmetry or vanishing first derivative, ‘a’ implies antisymmetry or vanishing value, and ‘e2’ implies extrapolated value.

The independent variables solved for the 1D accretion disc models presented in figures 7 and 8 are, $u_{x}$, $u_{y}$, $u_{z}$, $\ln\rho$, $s$ and $\nu_{\rm shock}$. The boundary conditions supplied for these 6 variables along the $z$-direction is given by bcz=‘s’,‘s’,‘a’,‘s:a2’,‘s:a2’,‘s’; where ‘a2’ implies antisymmetry with a vanishing second derivative.

The independent variables solved for the 3D accretion disc model presented in figures 9–11 are, $u_{x}$, $u_{y}$, $u_{z}$, $\ln\rho$, $s$ and $\nu_{\rm shock}$. The boundary conditions supplied for these 6 variables along the $z$-direction is given by bcz=‘s’,‘s’,‘a’,‘a2’,‘a2’,‘s’, while periodic boundary conditions are used along $x$ and $y$.

Appendix D Relation between pressure and density scale heights

As stated in section 4.3, $H_{\!p}=H_{\!\rho}$ for an isothermal stratification. For an isentropic stratification, we have $\gamma H_{\!p}=H_{\!\rho}$. In terms of the double-logarithmic temperature gradient, $\nabla={\rm d}{}\ln T/{\rm d}{}\ln p$, the general relation is given by

[TABLE]

Thus, since $\nabla_{\rm ad}=1-1/\gamma$, we have

[TABLE]

which is independent of $\gamma$. In the absence of convection, and for simple power law opacities given by (28) with single exponents $a$ and $b$, $\nabla=1/(1+n)$ depends on the polytropic index $n=(3-b)/(1+a)$ (Barekat and Brandenburg, 2014). In table LABEL:ScaleHeightRatio, some examples are listed.