Astrophysical simulations of convection frequently impose different thermal boundary conditions at the top and the bottom of the domain in an effort to more accurately model natural systems. In this work, we study Rayleigh-Benard convection (RBC) under the Boussinesq approximation. We examine simulations with mixed temperature boundary conditions in which the flux is fixed at the bottom boundary and the temperature is fixed at the top ("FT"). We aim to understand how FT boundaries change the nature of the convective solution compared to the traditional choice of thermal boundaries, in which the temperature is fixed at the top and bottom of the domain ("TT"). We demonstrate that the timescale of thermal relaxation for FT simulations is dependent upon the initial conditions. "Classic" initial conditions which employ a hydrostatically -- and thermally -- balanced linear temperature profile exhibit a long thermal relaxation. This long relaxation is not seen in FT simulations which use a TT simulation's nonlinear state as initial conditions ("TT-to-FT"). In the thermally relaxed, statistically stationary state, the mean behavior of an FT simulation corresponds to an equivalent simulation with TT boundaries, and time- and volume-averaged flow statistics like the Nusselt number and the Peclet number are indistinguishable between FT and TT simulations. FT boundaries are fundamentally asymmetric, and we examine the asymmetries that these boundaries produce in the flow. We find that the fixed-flux boundary produces more extreme temperature events than the fixed-temperature boundary. However, these near-boundary asymmetries do not measurably break the symmetry in the convective interior. We briefly explore rotating RBC to demonstrate that our findings with respect to thermal relaxation carry over to this more complex case, and to show the power of TT-to-FT initial conditions.

Numerical solutions of partial differential equations enable a broad range of scientific research. The Dedalus Project is a flexible, open-source, parallelized computational framework for solving general partial differential equations using spectral methods. Dedalus translates plain-text strings describing partial differential equations into efficient solvers. This paper details the numerical method that enables this translation, describes the design and implementation of the codebase, and illustrates its capabilities with a variety of example problems. The numerical method is a first-order generalized tau formulation that discretizes equations into banded matrices. This method is implemented with an object-oriented design. Classes for spectral bases and domains manage the discretization and automatic parallel distribution of variables. Discretized fields and mathematical operators are symbolically manipulated with a basic computer algebra system. Initial value, boundary value, and eigenvalue problems are efficiently solved using high-performance linear algebra, transform, and parallel communication libraries. Custom analysis outputs can also be specified in plain text and stored in self-describing portable formats. The performance of the code is evaluated with a parallel scaling benchmark and a comparison to a finite-volume code. The features and flexibility of the codebase are illustrated by solving several examples: the nonlinear Schrodinger equation on a graph, a supersonic magnetohydrodynamic vortex, quasigeostrophic flow, Stokes flow in a cylindrical annulus, normal modes of a radiative atmosphere, and diamagnetic levitation. The Dedalus code and the example problems are available online at http://dedalus-project.org/.

The magnetorotational instability (MRI) occurs when a weak magnetic field destabilises a rotating, electrically conducting fluid with inwardly increasing angular velocity. The MRI is essential to astrophysical disk theory where the shear is typically Keplerian. Internal shear layers in stars may also be MRI unstable, and they take a wide range of profiles, including near-critical. We show that the fastest growing modes of an ideal magnetofluid are three-dimensional provided the shear rate, $S$, is near the two-dimensional onset value, $S_c$. For a Keplerian shear, three-dimensional modes are unstable above $S\approx0.10S_c$, and dominate the two-dimensional modes until $S\approx2.05S_{c}$. These three-dimensional modes dominate for shear profiles relevant to stars and at magnetic Prandtl numbers relevant to liquid-metal laboratory experiments. Significant numbers of rapidly growing three-dimensional modes remain well past $2.05S_{c}$. These finding are significant in three ways. First, weakly nonlinear theory suggests that the MRI saturates by pushing the shear rate to its critical value. This can happen for systems, like stars and laboratory experiments, that can rearrange their angular velocity profiles. Second, the non-normal character and large transient growth of MRI modes should be important whenever three-dimensionality exists. Finally, three-dimensional growth suggests direct dynamo action driven from the linear instability.

A new gridding technique for the solution of partial differential equations in cubical geometry is presented. The method is based on volume penalization, allowing for the imposition of a cubical geometry inside of its circumscribing sphere. By choosing to embed the cube inside of the sphere, one obtains a discretization that is free of any sharp edges or corners. Taking full advantage of the simple geometry of the sphere, spectral bases based on spin-weighted spherical harmonics and Jacobi polynomials, which properly capture the regularity of scalar, vector and tensor components in spherical coordinates, can be applied to obtain moderately efficient and accurate numerical solutions of partial differential equations in the cube. This technique demonstrates the advantages of these bases over other methods for solving PDEs in spherical coordinates. We present results for a test case of incompressible hydrodynamics in cubical geometry: Rayleigh-B\'enard convection with fully Dirichlet boundary conditions. Analysis of the simulations provides what is, to our knowledge, the first result on the scaling of the heat flux with the thermal forcing for this type of convection in a cube in a sphere.

The Rossby number is a crucial parameter describing the degree of rotational constraint on the convective dynamics in stars and planets. However, it is not an input to computational models of convection but must be measured ex post facto. Here, we report the discovery of a new quantity, the Predictive Rossby number, which is both tightly correlated with the Rossby number and specified in terms of common inputs to numerical models. The Predictive Rossby number can be specified independent of Rayleigh number, allowing suites of numerical solutions to separate the degree of rotational constraint from the strength of the driving of convection. We examine the scaling of convective transport in terms of the Nusselt number and the degree of turbulence in terms of the Reynolds number of the flow, and we find scaling laws nearly identical to those in nonrotational convection at low Rossby number when the Predictive Rossby number is held constant. Finally, we describe the boundary layers as a function of increasing turbulence at constant Rossby number.

This paper demonstrates an equivalence between rotating magnetised shear flows and a stressed elastic beam. This results from finding the same form of dynamical equations after an asymptotic reduction of the axis-symmetric magnetorotational instability (MRI) under the assumption of almost-critical driving. The analysis considers the MRI dynamics in a non-dissipative near-equilibrium regime. Both the magnetic and elastic systems reduce to a simple one-dimensional wave equation with a nonlocal nonlinear feedback. Under transformation, the equation comprises a large number of mean-field interacting Duffing oscillators. This system was the first proven example of a strange attractor in a partial differential equation. Finding the same reduced equation in two natural applications suggests the model might result from other applications and could fall into a universal class based on symmetry.

This paper presents a method for the accurate and efficient computations on scalar, vector and tensor fields in three-dimensional spherical polar coordinates. The methods uses spin-weighted spherical harmonics in the angular directions and rescaled Jacobi polynomials in the radial direction. For the 2-sphere, spin-weighted harmonics allow for automating calculations in a fashion as similar to Fourier series as possible. Derivative operators act as wavenumber multiplication on a set of spectral coefficients. After transforming the angular directions, a set of orthogonal tensor rotations put the radially dependent spectral coefficients into individual spaces each obeying a particular regularity condition at the origin. These regularity spaces have remarkably simple properties under standard vector-calculus operations, such as \textit{grad} and \textit{div}. We use a hierarchy of rescaled Jacobi polynomials for a basis on these regularity spaces. It is possible to select the Jacobi-polynomial parameters such that all relevant operators act in a minimally banded way. Altogether, the geometric structure allows for the accurate and efficient solution of general partial differential equations in the unit ball.

We present a simulation code which can solve broad ranges of partial differential equations in a full sphere. The code expands tensorial variables in a spectral series of spin-weighted spherical harmonics in the angular directions and a scaled Jacobi polynomial basis in the radial direction, as described in Part-I. Nonlinear terms are calculated by transforming from the coefficients in the spectral series to the value of each quantity on the physical grid, where it is easy to calculate products and perform other local operations. The expansion makes it straightforward to solve equations in tensor form (i.e., without decomposition into scalars). We propose and study several unit tests which demonstrate the code can accurately solve linear problems, implement boundary conditions, and transform between spectral and physical space. We then run a series of benchmark problems proposed in Marti et al (2014), implementing the hydrodynamic and magnetohydrodynamic equations. We are able to calculate more accurate solutions than reported in Marti et al 2014 by running at higher spatial resolution and using a higher-order timestepping scheme. We find the rotating convection and convective dynamo benchmark problems depend sensitively on details of timestepping and data analysis. We also demonstrate that in low resolution simulations of the dynamo problem, small changes in a numerical scheme can lead to large changes in the solution. To aid future comparison to these benchmarks, we include the source code used to generate the data, as well as the data and analysis scripts used to generate the figures.

As the Science Mission Directorate contemplates establishing an open code policy, we consider it timely to share our experiences as the developers of the open-source partial differential equation solver Dedalus. Dedalus is a flexible framework for solving partial differential equations. Its development team primarily uses it for studying stellar and planetary astrophysics. Dedalus was developed originally for astrophysical fluid dynamics (AFD), though it has found a much broader user base, including applied mathematicians, plasma physicists, and oceanographers. Here, we will focus on issues related to open-source software from the perspective of AFD. We use the term AFD with the understanding that astrophysics simulations are inherently multi-physics: fluid dynamics coupled with some combination of gravitational dynamics, radiation transfer, relativity, and magnetic fields. In practice, a few well-known open-source simulation packages represent a large fraction of published work in the field. However, we will argue that an open-code policy should encompass not just these large simulation codes, but also the input files and analysis scripts. It is our interest that NASA adopt an open-code policy because without it, reproducibility in computational science is needlessly hampered.

Asteroseismology probes the interiors of stars by studying oscillation modes at a star's surface. Although pulsation spectra are well understood for solar-like oscillators, a substantial fraction of red giant stars observed by Kepler exhibit abnormally low-amplitude dipole oscillation modes. Fuller et al. (2015) suggests this effect is produced by strong core magnetic fields that scatter dipole internal gravity waves (IGWs) into higher multipole IGWs or magnetic waves. In this paper, we study the interaction of IGWs with a magnetic field to test this mechanism. We consider two background stellar structures: one with a uniform magnetic field, and another with a magnetic field that varies both horizontally and vertically. We derive analytic solutions to the wave propagation problem and validate them with numerical simulations. In both cases, we find perfect conversion from IGWs into magnetic waves when the IGWs propagate into a region exceeding a critical magnetic field strength. Downward propagating IGWs cannot reflect into upward propagating IGWs because their vertical wavenumber never approaches zero. Instead, they are converted into upward propagating slow (Alfvenic) waves, and we show they will likely dissipate as they propagate back into weakly magnetized regions. Therefore, strong internal magnetic fields can produce dipole mode suppression in red giants, and gravity modes will likely be totally absent from the pulsation spectra of sufficiently magnetized stars.

The annihilation of TeV photons from extragalactic TeV sources and the extragalactic background light produces ultrarelativistic $e^{\pm}$ beams, which are subject to powerful plasma instabilities that sap their kinetic energy. Here we study the linear phase of the plasma instabilities that these pair beams drive. To this end, we calculate the linear growth rate of the beam plasma and oblique instability in the electrostatic approximation in both the reactive and kinetic regimes, assuming a Maxwell-J{\"u}ttner distribution for the pair beam. We reproduce the well-known reactive and kinetic growth rates for both the beam plasma and oblique mode. We demonstrate for the oblique instability that there is a broad spectrum of unstable modes that grow at the maximum rate for a wide range of beam temperatures and wave vector orientations relative to the beam. We also delineate the conditions for applicability for the reactive and kinetic regimes and find that the beam plasma mode transitions to the reactive regime at a lower Lorentz factor than the oblique mode due to a combination of their different scalings and the anisotropy of the velocity dispersions. Applying these results to the ultrarelativistic $e^{\pm}$ beams from TeV blazars, we confirm that these beams are unstable to both, the kinetic oblique mode and the reactive beam-plasma mode. These results are important in understanding how powerful plasma instabilities may sap the energy of the ultrarelativistic $e^{\pm}$ beams as they propagate through intergalactic space.

It has been proposed that mixing induced by convective overshoot can disrupt the inward propagation of carbon deflagrations in super-asymptotic giant branch stars. To test this theory, we study an idealized model of convectively bounded carbon flames with 3D hydrodynamic simulations of the Boussinesq equations using the pseudospectral code Dedalus. Because the flame propagation timescale is much longer than the convection timescale, we approximate the flame as fixed in space, and only consider its effects on the buoyancy of the fluid. By evolving a passive scalar field, we derive a {\it turbulent} chemical diffusivity produced by the convection as a function of height, $D_{\rm t}(z)$. Convection can stall a flame if the chemical mixing timescale, set by the turbulent chemical diffusivity, $D_{\rm t}$, is shorter than the flame propagation timescale, set by the thermal diffusivity, $\kappa$, i.e., when $D_{\rm t}>\kappa$. However, we find $D_{\rm t}<\kappa$ for most of the flame because convective plumes are not dense enough to penetrate into the flame. Extrapolating to realistic stellar conditions, this implies that convective mixing cannot stall a carbon flame and that "hybrid carbon-oxygen-neon" white dwarfs are not a typical product of stellar evolution.

Spectral methods are an efficient way to solve partial differential equations on domains possessing certain symmetries. The utility of a method depends strongly on the choice of spectral basis. In this paper we describe a set of bases built out of Jacobi polynomials, and associated operators for solving scalar, vector, and tensor partial differential equations in polar coordinates on a unit disk. By construction, the bases satisfy regularity conditions at r=0 for any tensorial field. The coordinate singularity in a disk is a prototypical case for many coordinate singularities. The work presented here extends to other geometries. The operators represent covariant derivatives, multiplication by azimuthally symmetric functions, and the tensorial relationship between fields. These arise naturally from relations between classical orthogonal polynomials, and form a Heisenberg algebra. Other past work uses more specific polynomial bases for solving equations in polar coordinates. The main innovation in this paper is to use a larger set of possible bases to achieve maximum bandedness of linear operations. We provide a series of applications of the methods, illustrating their ease-of-use and accuracy.

The nonlinear evolution of the Kelvin-Helmholtz instability is a popular test for code verification. To date, most Kelvin-Helmholtz problems discussed in the literature are ill-posed: they do not converge to any single solution with increasing resolution. This precludes comparisons among different codes and severely limits the utility of the Kelvin-Helmholtz instability as a test problem. The lack of a reference solution has led various authors to assert the accuracy of their simulations based on ad-hoc proxies, e.g., the existence of small-scale structures. This paper proposes well-posed Kelvin-Helmholtz problems with smooth initial conditions and explicit diffusion. We show that in many cases numerical errors/noise can seed spurious small-scale structure in Kelvin-Helmholtz problems. We demonstrate convergence to a reference solution using both Athena, a Godunov code, and Dedalus, a pseudo-spectral code. Problems with constant initial density throughout the domain are relatively straightforward for both codes. However, problems with an initial density jump (which are the norm in astrophysical systems) exhibit rich behavior and are more computationally challenging. In the latter case, Athena simulations are prone to an instability of the inner rolled-up vortex; this instability is seeded by grid-scale errors introduced by the algorithm, and disappears as resolution increases. Both Athena and Dedalus exhibit late-time chaos. Inviscid simulations are riddled with extremely vigorous secondary instabilities which induce more mixing than simulations with explicit diffusion. Our results highlight the importance of running well-posed test problems with demonstrated convergence to a reference solution. To facilitate future comparisons, we include the resolved, converged solutions to the Kelvin-Helmholtz problems in this paper in machine-readable form.

Water's density maximum at 4C makes it well suited to study internal gravity wave excitation by convection: an increasing temperature profile is unstable to convection below 4C, but stably stratified above 4C. We present numerical simulations of a water-like fluid near its density maximum in a two dimensional domain. We successfully model the damping of waves in the simulations using linear theory, provided we do not take the weak damping limit typically used in the literature. In order to isolate the physical mechanism exciting internal waves, we use the novel spectral code Dedalus to run several simplified model simulations of our more detailed simulation. We use data from the full simulation as source terms in two simplified models of internal wave excitation by convection: bulk excitation by convective Reynolds stresses, and interface forcing via the mechanical oscillator effect. We find excellent agreement between the waves generated in the full simulation and the simplified simulation implementing the bulk excitation mechanism. The interface forcing simulations over excite high frequency waves because they assume the excitation is by the "impulsive" penetration of plumes, which spreads energy to high frequencies. However, we find the real excitation is instead by the "sweeping" motion of plumes parallel to the interface. Our results imply that the bulk excitation mechanism is a very accurate heuristic for internal wave generation by convection.

Thermal conduction is an important energy transfer and damping mechanism in astrophysical flows. Fourier's law - the heat flux is proportional to the negative temperature gradient, leading to temperature diffusion - is a well-known empirical model of thermal conduction. However, entropy diffusion has emerged as an alternative thermal conduction model, despite not ensuring the monotonicity of entropy. This paper investigates the differences between temperature and entropy diffusion for both linear internal gravity waves and weakly nonlinear convection. In addition to simulating the two thermal conduction models with the fully compressible Navier-Stokes equations, we also study their effects in the reduced, "sound-proof" anelastic and pseudo-incompressible equations. We find that in the linear and weakly nonlinear regimes, temperature and entropy diffusion give quantitatively similar results, although there are some larger errors in the pseudo-incompressible equations with temperature diffusion due to inaccuracies in the equation of state. Extrapolating our weakly nonlinear results, we speculate that differences between temperature and entropy diffusion might become more important for strongly turbulent convection.

The speed of sound greatly exceeds typical flow velocities in many stellar and planetary interiors. To follow the slow evolution of subsonic motions, various sound-proof models attempt to remove fast acoustic waves whilst retaining stratified convection and buoyancy dynamics. In astrophysics, anelastic models typically receive the most attention in the class of sound-filtered stratified models. Generally, anelastic models remain valid in nearly adiabatically stratified regions like stellar convection zones, but may break down in strongly sub-adiabatic, stably stratified layers common in stellar radiative zones. However, studying stellar rotation, circulation, and dynamos requires understanding the complex coupling between convection and radiative zones, and this requires robust equations valid in both regimes. Here we extend the analysis of equation sets begun in Brown Vasil & Zweibel 2012, which studied anelastic models, to two types of pseudo-incompressible models. This class of models has received attention in atmospheric applications, and more recently in studies of white-dwarf supernovae progenitors. We demonstrate that one model conserves energy but the other does not. We use Lagrangian variational methods to extend the energy conserving model to a general equation of state, and dub the resulting equation set the Generalized Pseudo-Incompressible (GPI) model. We show that the GPI equations suitably capture low frequency phenomena in both convection and radiative zones in stars and other stratified systems, and we provide recommendations for converting low-Mach number codes to this equation set.

Typical flows in stellar interiors are much slower than the speed of sound. To follow the slow evolution of subsonic motions, various sound-proof equations are in wide use, particularly in stellar astrophysical fluid dynamics. These low-Mach number equations include the anelastic equations. Generally, these equations are valid in nearly adiabatically stratified regions like stellar convection zones, but may not be valid in the sub-adiabatic, stably stratified stellar radiative interiors. Understanding the coupling between the convection zone and the radiative interior is a problem of crucial interest and may have strong implications for solar and stellar dynamo theories as the interface between the two, called the tachocline in the Sun, plays a crucial role in many solar dynamo theories. Here we study the properties of gravity waves in stably-stratified atmospheres. In particular, we explore how gravity waves are handled in various sound-proof equations. We find that some anelastic treatments fail to conserve energy in stably-stratified atmospheres, instead conserving pseudo-energies that depend on the stratification, and we demonstrate this numerically. One anelastic equation set does conserve energy in all atmospheres and we provide recommendations for converting low-Mach number anelastic codes to this set of equations.

Recently, Silvers, Vasil, Brummell, & Proctor (2009), using numerical simulations, confirmed the existence of a double diffusive magnetic buoyancy instability of a layer of horizontal magnetic field produced by the interaction of a shear velocity field with a weak vertical field. Here, we demonstrate the longer term nonlinear evolution of such an instability in the simulations. We find that a quasi two-dimensional interchange instability rides (or "surfs") on the growing shear-induced background downstream field gradients. The region of activity expands since three-dimensional perturbations remain unstable in the wake of this upward-moving activity front, and so the three-dimensional nature becomes more noticeable with time.

In this paper we discuss the effect of recombinations to highly excited states (n > 100) in hydrogen during the cosmological recombination epoch. For this purpose, we developed a new ODE solver for the recombination problem, based on an implicit Gear's method. This solver allows us to include up to 350 l-resolved shells or ~61 000 separate levels in the hydrogen model and to solve the recombination problem for one cosmology in ~27 hours. This is a huge improvement in performance over our previous recombination code, for which a 100-shell computation (5050 separate states) already required ~150 hours on a single processor. We show that for 350 shells down to redshift z ~200 the results for the free electron fraction have practically converged. The final modification in the free electron fraction at z ~200 decreases from about \DeltaNe/Ne ~2.8% for 100 shells to \DeltaNe/Ne ~1.6% for 350 shells. However, the associated changes in the CMB power spectra at large multipoles l are rather small, so that for accurate computations in connection with the analysis of Planck data already ~100 shells are expected to be sufficient. Nevertheless, the total value of \tau could still be affected at a significant level. We also briefly investigate the effect of collisions on the recombination dynamics. With our current estimates for the collisional rates we find a correction of \DeltaNe/Ne ~ -0.088% at z ~ 700, which is mainly caused by l-changing collisions with protons. Furthermore, we present results on the cosmological recombination spectrum, showing that at low frequencies collisional processes are important. However, the current accuracy of collisional rates is insufficient for precise computations of templates for the recombination spectrum at \nu<~1 GHz, and also the effect of collisions on the recombination dynamics suffers from the uncertainty in these rates.

Previous theoretical work has speculated about the existence of double-diffusive magnetic buoyancy instabilities of a dynamically evolving horizontal magnetic layer generated by the interaction of forced vertically sheared velocity and a background vertical magnetic field. Here we confirm numerically that if the ratio of the magnetic to thermal diffusivities is sufficiently low then such instabilities can indeed exist, even for high Richardson number shear flows. Magnetic buoyancy may therefore occur via this mechanism for parameters that are likely to be relevant to the solar tachocline, where regular magnetic buoyancy instabilities are unlikely.

The turbulent diffusion tensor describing the evolution of the mean concentration of a passive scalar is investigated for non-helically forced turbulence in the presence of rotation or a magnetic field. With rotation, the Coriolis force causes a sideways deflection of the flux of mean concentration. Within the magnetohydrodynamics approximation there is no analogous effect from the magnetic field because the effects on the flow do not depend on the sign of the field. rotation and magnetic fields tend to suppress turbulent transport, but this suppression is weaker in the direction along the magnetic field. Turbulent transport along the rotation axis is not strongly affected by rotation, except on shorter length scales, i.e., when the scale of the variation of the mean field becomes comparable with the scale of the energy-carrying eddies. These results are discussed in the context of anisotropic convective energy transport in the Sun.