© ESO, 2016

1. Introduction

Stellar evolution modeling is subject to various uncertainties related to (magneto)hydrodynamic processes, such as mixing of chemical elements and the transport of angular momentum. In 1D models, these uncertainties are represented by a parameter that is then usually calibrated against observations. However, values calibrated from one star do not always fit other stars, and there can be degeneracy when several parameters are used. A good example is convection, which is normally parametrized in terms of a mixing length and an overshooting parameter. We cannot predict the appropriate values from theory, presumably because mixing length theory (MLT) is a physically incorrect model: experiments and simulations of convection do not show eddies mixing over some lengths, but rather upward- and downward-moving fluid packets surviving over large distances, whereby the boundaries that are absent from MLT are critical (see e.g., Stein & Nordlund 2006). Unfortunately though, there is not yet any convenient but physically more realistic way to parametrize convection in a 1D model. Other uncertainties come from our lack of understanding of the physics of stellar winds and mass loss, and in binaries: the physics of accretion, mergers, tides and magnetic braking.

The numerical scheme described here is designed for the investigation of physical processes in stellar radiative zones. It is known that both chemical elements and angular momentum can be transported through radiative zones, but the physics is poorly understood. Mixing must result from fluid motions in the radial direction, despite the fact that radial motion is suppressed in a radiative zone by the strong buoyant restoring force. Magnetic torques probably play a crucial role, especially in the transport of angular momentum. In any case, flows are expected to be predominantly horizontal, with a relatively small radial component. The numerical scheme described here is specifically designed for these predominantly horizontal flows with an extreme horizontal-vertical aspect ratio. In the next section we describe some physical processes that this code can be used to examine.

2. Physical processes in radiative zones

Many interesting effects in radiative zones come from rotation. Rotation makes the star oblate, and this inevitably means that the gradients of thermodynamical variables and gravitational potential are no longer parallel (Von Zeipel 1924). Eddington (1925), Vogt (1925), and Sweet (1950) suggested that this gives rise to a meridional flow, but then Busse (1981) pointed out that while such a flow satisfies the heat equation, it does not satisfy the momentum equation. A physical solution was described by Zahn (1992) in terms of the baroclinic instability and the predominantly horizontal turbulence in which it results. This gave rise to the concept of “shellular rotation”, where latitudinal variations in angular velocity Ω are quickly dissipated, leaving only radial variation and Ω = Ω(r). In addition, latitudinal variations in composition are smoothed. This concept, while based on sound physical reasoning, was welcomed by the stellar evolution community also because of its convenience: it meant that work could continue with 1D modeling. Although nobody has found any holes in this elegant theory, it seems to produce insufficient transport of angular momentum for reproducing for instance, the slow rotation of white dwarfs (e.g., Suijs et al. 2008; see also Meynet & Maeder 1997). There may be some missing ingredient, or perhaps some modification is required; in any case it is certainly a useful exercise to test the physical model with simulations.

In differentially rotating radiative zones, one also expects magnetic fields to be important. Differential rotation will wind up any seed field into a predominantly toroidal configuration, which is then unstable to an interchange instability (Tayler 1973), and in so decaying, a new poloidal field is produced that can be wound up anew by the differential rotation, thereby closing a dynamo loop (Spruit 2002). This dynamo process should result in transport of angular momentum and chemical elements in the radial direction, and a prescription for this has been implemented in some stellar evolution models (e.g., Heger et al. 2005; Suijs et al. 2008). However, Spruit’s theory suffers from the same problem as Zahn’s: although based on sound physical reasoning, it has never been tested rigorously. It is popular mainly on account of its convenience for implementation in 1D stellar evolution models (as a “magic formula”). Although it is currently the only model that produces anywhere near the angular momentum transport necessary to explain, for instance, the slow rotation of white dwarfs, there are quantitative problems reproducing, for instance, the core rotation of red giant stars as measured via asteroseismology with the Kepler satellite (Mosser et al. 2012). Finally, there have been claims that the Tayler-Spruit mechanism cannot work because of the lack of an axisymmetric mode; it turns out, however, that this is both misleading and untrue, because (a) there is indeed a strong axisymmetric mode and (b) an axisymmetric mode is not necessary anyway. See Ibáñez-Mejía & Braithwaite (2015) for a discussion of this point. Braithwaite (2006) and Rüdiger et al. (2015) find the mechanism in simulations. However, further numerical work is needed.

Finally, note that transport of angular momentum in radiative zones may also have to do with the asymmetric dissipation of internal gravity waves (Talon & Charbonnel 2005; Barker & Ogilvie 2010; Alvan et al. 2013). The numerical scheme described here is unfortunately unsuitable for investigating this, although this subtle effect seems so far to have been difficult to find with other schemes.

3. Numerical modeling

In hydrodynamical (HD) and magnetohydrodynamical (MHD) simulations, a set of partial differential equations is numerically integrated forward in time, in small increments called the “time step”. The time step is subject to a range of limits to do with the speed of propagation of information; numerical schemes generally become unstable if information is allowed to propagate further than some fraction of a grid spacing in one time step. For instance, a simple Cartesian hydrodynamical code might have the following time step: (1)where Δx, Δy and Δz are the grid spacings in the three dimensions, u is the gas velocity, c_s is the sound speed and C is some dimensionless constant whose value will depend on properties of the numerical discretisation scheme, but might be for instance 0.5. This would ensure that no information can propagate more than 0.5 grid spacings in any direction during one time step. Other codes will have additional, similar restrictions from the propagation of Alfvén waves and from diffusion; for instance, the sound speed c_s in the expression above might be replaced by the fast magnetosonic speed.

It is often possible to make approximations that remove some modes of information propagation, allowing a larger time step. In addition to the basic constant-density incompressible approximation, in which sound and buoyancy waves are both absent, the anelastic (Ogura & Phillips 1962) and Boussinesq approximations are commonplace (see e.g., Lilly 1996 for a review and comparison). Both are widely used in astrophysical hydrodynamics, for instance in studies of convection in planetary and stellar interiors (e.g., Browning 2008; Chen & Glatzmaier 2005).

In gravitationally-stratified fluids in which the characteristic length scale of the flow is much greater in the horizontal direction than in the vertical, modeling requires Δz ≪ Δx, Δy. In these contexts it can be fruitful to make the approximation of a hydrostatic equilibrium, in which perfect vertical force balance is assumed (2)In this approximation, the fluid adjusts to vertical force balance on the timescale of sound waves propagating in the vertical direction Richardson (1922), and the approximation is valid as long as this is much shorter than all other timescales of interest. Vertically propagating sound waves are filtered out (as well as high frequency internal gravity waves) and the z-component of the time step restriction (Eq. (1)) can be removed. For gravitationally-stratified flows with an extreme aspect ratio such as those expected in a radiative zone, this approximation makes more sense than the anelastic or Boussinesq approximations since it is “closer to the truth” in the sense of representing a smaller compromise in terms of accuracy. It also has the advantage over these two approximations that one is not restricted to modeling gentle phenomena close to a constant background state. It allows, for instance, supersonic (in the horizontal direction), compressible flow.

The hydrostatic approximation reduces the number of independent variables in the system. For instance, with “raw” hydrodynamic equation, there are two thermodynamic variables plus three components of velocity. In the equivalent hydrostatic system, the vertical component of the velocity is no longer independent, but calculated by integration of Eq. (2); furthermore one of the thermodynamic variables is lost, as is described in the next section.

In the astrophysical context we often want to model conducting fluids with magnetic fields. Note that although the hydrostatic approximation filters out vertically propagating sound waves, vertically propagating magnetic waves are not entirely filtered; a magnetohydrostatic scheme is therefore of use only in the case where the plasma β is high, i.e., where the Alfvén speed is much less than the sound speed. In practice this is usually not an additional restriction, since essentially all astrophysical contexts where the hydrostatic approximation is useful (i.e., fulfill the aforementioned conditions) have high β.

There are various ways in which the hydrostatic approximation can be implemented, resulting in different sets of equations and independent variables. The vertical coordinate can be height, pressure, entropy or some combination of these, and the upper and lower boundaries can be fixed in height, pressure, entropy or some combination (see e.g., Kasahara 1974 or Konor & Arakawa 1997, for a review of coordinate systems). Not surprisingly, the choices of coordinates and boundaries are not independent of each other, but there are a number of possible combinations. The best combination depends on the context. For instance, Braithwaite & Cavecchi (2012) describe a numerical scheme using the pressure hybrid σ coordinate, with lower and upper boundaries fixed in height and pressure respectively. This scheme is particularly suitable for neutron star oceans because the upper boundary can move upwards, allowing for large temperature variations, and has recently been used to model Type I X-ray bursts (Cavecchi et al. 2013, 2015).

Kasahara & Washington (1967) describe a hydrostatic numerical scheme using height coordinates which they used to construct a global circulation model of the atmosphere, in one of the earliest examples of 3D numerical hydrodynamics. However, owing to advances in hybrid coordinate systems which allowed also for topographical features – mountain ranges and so on – this scheme enjoyed only brief popularity in atmospheric physics. In the context of stellar interiors, however, the requirements are different. To start, the ideal scheme will have the same type of upper and lower boundary. There are two obvious choices: height coordinates, with boundaries fixed in height, and pressure coordinates, with boundaries fixed in pressure. Both have their pros and cons. In pressure coordinates, the equations are a somewhat simpler. In hydrostatic equilibrium, the pressure at any point is simply equal to the weight of the column of gas above that point and with pressure coordinates each grid box (which has a constant pressure difference ΔP from top to bottom) will contain constant mass. The continuity equation therefore becomes (3)which has the same form as the familiar incompressible equation ∇·u = 0 where vertical velocity u_z ≡ dz/ dt has been replaced by ω ≡ dP/ dt. One could also, in radiative zones, consider entropy coordinates; these have the advantage that the vertical “velocity” ds/ dt is small, a function only of heating and cooling, which reduces numerical diffusion in the vertical direction. However, although the system with the height coordinate has equations which look more complex, this coordinate has two significant advantages. The first is simply that the implementation of magnetic fields is more complicated in both pressure and entropy coordinates than in height coordinates: it is obviously more complicated to calculate for instance the curl of a vector when the three coordinates dont have the same units. Secondly, the equations only need to be integrated in the vertical direction, whereas with pressure coordinates it is necessary to solve the Poisson equation (∇²a = b where a is the unknown). This is not as straightforward as it looks, and is tricky to parallelize efficiently. (Note that this is also a disadvantage with Boussinesq and anelastic schemes.) It is largely for these reasons that we choose height coordinates to investigate stellar interiors.

In Sect. 4 we present the basic equations, presenting simple test cases in Sect. 5 and summarizing in Sect. 6.

4. Method

Our method is based on the equations developed by Kasahara & Washington (1967), which utilize the hydrostatic approximation, for use in atmospheric circulation modeling. We present here, a similar scheme which uses spherical coordinates and allows for the inclusion of magnetic fields, making it well suited for the modeling of various phenomena in gravitationally stratified fluids. In this section we describe the basic equations, coordinate system conventions and respective variables, the time stepping, the numerical grid and boundary conditions.

4.1. Basic equations

We use a similar symbol convention and scheme as Kasahara & Washington (1967), who also use spherical coordinates with height as the vertical (radial) coordinate. The symbols to be used are defined in Table 1.

In this set of coordinates, assuming hydrostatic equilibrium and a thin shell, the momentum conservation equation becomes (4)Here F_L,λ and F_L,ϕ are the horizontal components of the Lorentz force, and F_V,λ and F_V,ϕ are the components of the viscous force per unit volume. They are described below.

A discussion of the approximations used to reach this form is presented by Phillips (1966). Now, let us define the Lagrangian derivative, (5)The mass continuity equation is (6)We also assume an ideal gas equation of state, that is (7)The thermodynamic equation may be defined, using Eq. (5). (8)where a detailed explanation of the heating per unit mass Q is given in Sect. 4.2.

The hydrostatic approximation is imposed by assuming that there is hydrostatic equilibrium, i.e., Eq. (2). (9)Here, we are assuming that gravity is constant which is valid if the domain being simulated is a very thin shell. Integrating the thermodynamic equation, Eq. (8), from a vertical point z to the top of the simulation which we call z_T, where z_T is the radial domain in the simulation, we obtain (10)where at r = z_T + a so that the vertical speed u_z is 0 at r = z_T + a.

Now, using Eqs. (9), (10) and (5), one can write the Lagrangian derivative for the pressure: (11)where (12)Now if we put Eq. (11) into Eq. (8) and integrate between the bottom of the vertical domain (z = 0) to a point z, we can obtain the vertical speed u_r. Here we assume that u_r = 0 at z = 0, i.e., there is no flow across the lower boundary. (13)Now we can compute B so that the boundary condition for Eq. (13) is satisfied: (14)The next step is to define the Lorentz force (per unit volume): (15)The induction equation we use is (16)where η is the Ohmic diffusivity. In a numerical scheme where the grid boxes have extreme aspect ratios, it is not obvious what is the best way to calculate the diffusive term. This brings us onto the next subsection.

4.2. Diffusive scheme

Now, the viscous force terms in Eq. (4) may be expressed as (17)where ν represents the kinematic viscosity, k is a dimensionless scaling factor described below, and (18)is the horizontal Laplacian. We ignore terms proportional to the velocities and first derivatives, since they are small compared with the second derivatives.

We use a scaling factor k to make viscous time scales similar in the vertical and horizontal directions, (19)where Δr is the radial grid spacing, Δs is the minimum of the two horizontal grid spacings and K_v is an arbitrary constant of order of unity to disentangle the two viscous timescales, if required. Viscosity ν will be discussed later with the numeric scheme. Although there is no real physical justification for the anisotropy in the kinematic viscosity, it is convenient from a computational point of view, which will be more clear in the next section.

When the diffusive current is calculated in the induction equation, Eq. (16), the scale factor k from Eq. (19) is again used wherever a vertical derivative is present. This is done to ensure that the magnetic diffusion timescale in the much shorter radial direction is scaled to be similar to the longer horizontal directions, resulting in similar magnetic-diffusive time step restrictions.

The heating rate per unit mass Q in Eq. (8) and subsequent equations can be expressed as (20)Where the heating rate has been split into its constituent parts: heating due to thermal diffusion, viscous heating and Ohmic heating. The first can be expressed as (21)where κ represents the thermal diffusivity, k is the same scaling factor as with the kinematic viscosity and T is the temperature.

We use the standard viscous dissipation function φ as Q_visc: (22)where ∇u is the gradient of the velocity field, hence a tensor, u, the colon is the contraction of two tensors, and is the viscous stress tensor, (23)Here is the identity matrix, λ is the bulk viscosity coefficient, and is the strain rate tensor, defined as (24)We assume the bulk viscosity to be negligible. Using the same approximations to this set of coordinates to consider radial terms and the anisotropic vertical viscosity factor k, the heating rate due viscosity is (25)The Joule heating is computed as (26)where the same scale factor k is used with vertical derivatives.

4.3. Time stepping

We use a third order, 2-N storage Runge-Kutta scheme to advance our predicted variables in time. Several Courant conditions are imposed to limit information from moving more than a certain fraction of a grid spacing during one time step. From each of these conditions a time step is calculated; the minimum of these time steps is then used. The first condition concerns waves, which carry information at their speed relative to the fluid, plus the flow speed. In our scheme, this works differently in the horizontal and vertical. The resulting time steps are: where as before Δs is the minimum of the two horizontal grid spacings, and is the sound speed in an ideal gas, is the Alfvén speed, and the subscripts h and z refer to horizontal and vertical components of the variables. The difference between the expressions above comes from the fact that we have filtered out vertically propagating sound waves, but that Alfvén waves can still propagate vertically; in the horizontal direction the maximum wave speed is the magnetosonic speed. The dimensionless constant c₁ is generally set to less than 0.5; the optimum value is difficult to calculate from first principles and so is essentially a matter of trial and error.

Diffusion also transmits information. From dimensional analysis we find a time step (29)where ν_i is any of the three given diffusivities (momentum ν, thermal κ and magnetic η: whichever gives the lowest time step is used), and c₂ is another dimensionless constant similar to c₁ described above. Here, it is obvious that if we had similar diffusivities in the vertical and horizontal directions, then the much smaller Δr would restrict the whole calculation. That’s why we scaled previously the vertical viscosity and therefore now we can neglect the vertical dimension.

Finally, we don’t want our predicted variable P to vary too much over one time step, therefore we have, using another dimensionless constant c₃: (30)This standard approach helps ensure that the pressure, which is an independent variable, is not allowed to change too quickly, which if left unchecked can cause numerical problems. A similar restriction on the temperature might make sense in the future, but for the moment we do not envisage using the code for situations where the temperature can vary quickly. The final time step will be the minimum of these time steps at every grid point.

4.4. Diffusivities

The kinematic viscosity is (31)where n₁ is a dimensionless constant, c_max is the maximum flow + sound speed, and Δs is the minimum distance in the horizontal directions. This expression helps to disentangle the tunable parameter nu₁ to the more physical c_max, removing the need to assign units to the viscosity parameter. The formalism is applicable if c_max does not vary much over the domain and over time. We test this on a shear instability simulation.

The thermal and magnetic diffusivities are related to the kinematic viscosity by dimensionless constants, so that (32)where n₂ is the inverse of the Prandtl number and n₃ is the inverse of the magnetic Prandtl number by definition. This of course presumes no significant temperature or current structures.

To save computational time, there is an option in the code to use the assumption that variations in c_max are small and to use a constant viscosity, so that Δt_diff is also constant and is computed just at the beginning of the simulation, rather than having to calculate it every time step from Eq. (29).

4.5. Numerical grid

We utilize a staggered grid, using the Marker-and-Cell (MAC) method developed by Harlow & Welch (1965). We define P, ρ and T at the center of each grid cell, while the momenta and magnetic fields are defined at the respective faces, as shown in Fig. 1. To solve the equations, interpolations, derivatives and integrals are needed. We use a set of six points in each of these operations to get a sixth order interpolation and a fifth order derivative. A set of five points is used at each point to get a fifth order integral. The grid is regular in each direction, the center is the geometrical center of each grid cell. Appendix A contains the functions and derivations which we use.

Fig. 1 Positions in each grid cell at which various quantities are defined. The longitudinal, latitudinal and radial components of velocity are denoted by u, v and w.

4.6. Boundary conditions

In the radial boundaries, symmetric boundary conditions are used for the horizontal components of the momenta, while all other variables are calculated with the aid of ghost zones. In the ghost zones, we calculate the pressure by means of a logarithmic extrapolation, while the components of the magnetic fields can be chosen as needed. For the horizontal boundaries, we utilize periodic boundary conditions.

5. Results

In this section we describe various tests.

5.1. Plane wave

The first test to apply is a simple sound wave. We begin with an isothermal set of initial conditions, i.e., an exponential profile in pressure in the radial direction, and a fluid at rest. To test a purely 1D case, we set both the latitude and longitude ranges to be small, 0.1 radians. We also set a minimum radius, a, of 1000 length units, and a radial range of 0.1 length units.

To simulate the propagation of a pressure profile, or wave packet, we add to the pressure profile an overpressure constant in radius: a Gaussian profile in longitude with a full-width- at-half-maximum of 7 grid cells, and normalized to have a peak of the order of 0.1% of the average pressure at the bottom of the domain. The results can be seen in Fig. 2. The overpressure creates a pair of waves which travel to the edges and come back through the periodic boundaries. The wave should be non-dispersive, therefore the density profile after one sound crossing time should be exactly the same. This can be seen in Fig. 3, where the difference between the pressure after one sound crossing time and the initial pressure is plotted against longitude. The fractional difference is of order 10^-6; given that the relative amplitude of the wave is 10^-3, we see that the code has preserved the waveform with an accuracy of 10^-3. We find that the sound crossing time for the wave is 0.79 of the theoretical sound crossing time.

Fig. 2 Contour plot of the pressure in the plane wave test at the center of the radial range. The pressure perturbation moves towards the edges, preserving shape, like a 1D. After one period, the wave is back to its initial position. This behavior is constant in radius.

Fig. 3 Fractional difference between the pressure after one period and the initial pressure, taken at a height in the center of the domain. We can see that there is no clear dispersion of the wave packet, as expected from an inviscid fluid.

If we simulate a larger range in latitude, then the effect of the curvature becomes apparent. In Fig. 4 we see what happens when we have a latitude range of 2.0 radians. The initial setup is otherwise the same as the previous problem. Here, we find that the sound crossing time for the wave is 0.78 of the theoretical sound crossing time.

Fig. 4 Contour plot of the pressure in the plane wave test with higher angular range, with the expected theoretical position of a wave traveling with the local sound speed as a black dashed line. The wave appears to move faster away from the equator because the grid longitudinal spacing is smaller by a factor cosϕ.

5.2. 2D wave

Using the same background initial conditions and resolution as the previous 1D case with a short range in latitude and longitude, we add this time a 2D overpressure, constant in radius: a 2D Gaussian profile with a full-width-at-half-maximum of 7 grid cells, and normalized to have a peak of the order of 1% of the average pressure. Here we also switched off the viscosity. We can see in Fig. 5 the evolution of the system before the wave reaches the horizontal boundaries.

5.3. Shear instability

To test this instability, we simulated a box with the same dimensions as that in the previous section, with no radial dependence other than the hydrostatic equilibrium, and a discontinuous flow in the longitudinal direction with the following profile for the longitudinal velocity u: (33)where c_s is the sound speed, ϕ is the latitude and ϕ_T is the size of the box in latitude. Since the box is per construction symmetrical about the equator, this profile simulates a zone about the equator rotating as a rigid solid in one direction, and two zones counter-rotating. The speed difference is set at 10% of the sound speed to limit ourselves to the subsonic regime. We use a value of nu₁ = 0.005, a Prandtl number of 1, and no magnetic fields on this simulation.

Fig. 5 Contour plot of the pressure in the 2D wave test. We can see how it expands as a circular wave.

We add a perturbation in the longitudinal velocity of the form 0.01 | u_o | cos(2πnλ/λ_T), where λ is the longitude, and λ_T is the size of the box in longitude. n here is equivalent to the wavenumber. We used 5 values of wavenumbers, from 1 to 5, to test the growth rate of the instability.

We use a resolution of 50 × 256 × 256 in order to resolve the shorter wavelengths, giving us a resolution of about 51 × 51 grid cells per wavelength in the n = 5 case.

In Fig. 6 we can see the initial conditions of longitudinal velocity in the top-left, and a later snapshot of the longitudinal velocity for a system with an n = 5 perturbation in the bottom-left. The sinusoidal profile present in the longitudinal velocity as a perturbation is propagated to the latitudinal velocity as can be seen in the bottom-right. On the top-right side we can see what happens with the material with the help of a passive tracer. This variable acts as a dye, allowing us to follow material during the simulation. The governing equation of the passive tracer is a simple conservation equation, with a modest level of diffusion, but it does not interact in any other way with the fluid. Its initial conditions for this experiment is 0.5 in one zone, and −0.5 in the other.

The snapshot, taken after 3.6 sound crossing times (τ_s), where for the longitudinal domain λ_T and average sound speed , shows the typical spirals of a shear mixing through the passive tracer, and the different velocities.

Fig. 6 Shear instability in a square box of 50 × 256 × 256 resolution with a perturbation of n = 5. Top-left: initial longitudinal velocity u. The other three frames are at t = 281, approximately 3.6 sound crossing times: top-right: passive tracer used to follow the locations of the 2 fluids; bottom-left: the longitudinal velocity u and bottom-right: latitudinal velocity v. These horizontal cuts are taken half-way up the computational domain.

Now, the amplitude of the instabilities should grow exponentially, in the linear regime, as e^σt. The growth rate σ can be expressed, for an incompressible flow and fluids of equal density (e.g., Choudhuri 1998): (34)where Δu is the velocity difference across the discontinuity and k is the wavenumber. In this case where ϕ_disc is taken at the latitude of the discontinuity, and k = 2πn/ [ λacosϕ_disc ].

The simulations give growth rates fairly close to this prediction. We can see the exponential growth of the amplitudes of the first five wavenumbers of an initially n = 1−5 perturbed model, from a Fourier analysis done on the latitudinal component of the velocity in Fig. 7. The n = 1 mode (solid line) is the one which grows slowest, as expected, with faster growths for n = 2, 3, 4, and 5. Note that the growth rate is computed up to only 3 sound crossing times, since the instability leaves linearity when the amplitude is comparable to the wavelength of the perturbation. The growth rates do not agree entirely with the prediction in Eq. (34), which is presumably for two reasons: at high wavelength (e.g., n = 1), the wavelength is comparable to the dimension of the fluid in the direction perpendicular to the flow, whereas the prediction is calculated assuming the fluid has infinite extent in that direction. At low wavelength, the instability is damped by diffusion.

Here we can say something about the variation in the fluid speed plus the sound speed, hence the viscosity. The total variation of c_max over time is about ~0.4%. The maximum absolute difference between 2 individual values of c_max is ~14%.

Fig. 7 Amplitude and growth rate of shear instability modes. Top: amplitude as a function of time, in units of the sound crossing time 0.5τs. Bottom: the growth rate measured after t = 0.6τs (determined by taking the slope of a linear fit of the logarithm of the perturbation amplitudes), are roughly proportional to the wave number n, and roughly agree with the predicted growth rate from Eq. (34).

In addition to testing the shear instability due to a latitudinal shear, we tested the code against a radial shear. The code assumes a hydrostatic equilibrium in the radial direction, testing a radial shear instability is important. The setup of this experiment is similar to the longitudinal shear, with the exception that the longitudinal velocity experiences a discontinuity in the radial direction as:(35)where a is again the minimum value of the radius in the simulation and z_T is the radial extent of the simulation. As before, a velocity perturbation is added to initiate the shear instability. A perturbation to the longitudinal velocity, u, is added as: u = u_o + 0.01 | u_o | cos(2πnλ/λ_T), where once again λ is the longitude, and λ_T is the size of the box in longitude. n here is equivalent to the wavenumber. The simulation was run with a resolution of 96 × 96 × 32, and given an initial perturbation to each of the wavenumbers n = 1–5. Figure 8 shows an equatorial cut of the passive tracer at the start of the simulation and after 7.5 sound crossing times. We see a strong contribution from the n = 5 wavenumber, as it is the wavenumber which grows most quickly in the linear regime, as can be seen in Fig. 9.

Fig. 8 Shear instability in a square box with the shearing occurring in the radial direction for a model with an initial n = 1–5 perturbation. Panel a) shows the passive scalar at the start of the simulation in an equatorial cut. Panel b) shows the passive tracer after 7.5 sound crossing times, here it is evident that the n=5 mode dominated in the early evolution.

Fig. 9 Growth rate of the n = 2–5 modes for the shear instability when the shearing boundary is in the radial direction, for the first 1.8 sound crossing times. The n = 5 mode is the solid blue line, the n = 4 mode is the dashed gold line, the n = 3 mode is the dotted black line, and the n = 2 mode is the dashed-dotted red line.

5.4. Rotation

Our code does support a rotating coordinate system, therefore a natural test would be to test the propagation of a 2D wave with a non-zero rotational velocity Ω. The relevant parameter here is the Rossby number R_o, or inertial to Coriolis force ratio, which is defined as (36)This number depends on the relevant length L, the relevant velocity U and the latitude ϕ. In our simulations the typical value ~10^-3, similar to cyclones, which represents a system dominated by the Coriolis force. We use a value of nu₁ = 0.02, a Prandtl number of 1, and no magnetic fields on this simulation.

We use a grid of 50 × 128 × 256, with a latitude range of 2 radians, a longitude range of 1 radian, and a radial extent of 0.1 radial units with a minimum radius of a = 1000 radial units, which is roughly equivalent to 60% of the sphere. We simulate two systems at different absolute latitudes, one with a Gaussian overpressure and the other with a Gaussian under-pressure. Both perturbations have initially the same radius. A plot of pressure versus latitude in the initial conditions can be seen in Fig. 10.

The results can be seen in Fig. 11. Two snapshots are shown, one after just 20 time steps, where the fluid is starting to move, and another after 100, where the fluid is already rotating. Plotted are streamlines on a horizontal surface halfway up the computational domain. It is apparent that both systems rotate in the same direction, and that their radius is proportional to 1 / sinϕ.

Fig. 10 Pressure normalized to the surrounding pressure versus latitude at the longitude of the perturbation. The southern hemisphere perturbation is a overpressure relatively far from the equator, and the northern perturbation is relatively close to the equator.

5.5. Alfvén waves

We simulate a simple transverse magnetic wave in the vertical direction. We used a 300 × 16 × 16 grid, with a latitude domain of 0.1 radians, a longitude domain of 0.1 radians, and a radial range of 1 with a minimum radius of a = 1000. We set a perturbation in the longitudinal velocity u = max(0,u₀(1−r/r_p)) where u₀ is the maximum velocity of the perturbation, and r_p the radius at which the perturbation disappears. We use a value of nu₁ = 0.02, a Prandtl number of 1, and a magnetic Prandtl number of 1.

The results are shown as a time series in Fig. 12. The comparison with the predicted Alfvén velocity in Fig. 13 is very good, except where the wave starts and bounces and the velocity is more difficult to compute, since the wave is compressed at the boundaries.

Fig. 11 Streamlines for two snapshots of two rotating over- and under-pressurized features rotating and balanced by the Coriolis force, taken in the middle of the radial domain. The left panel shows how the fluid starts to move away from the center of the southern perturbation, and towards the center of the northern perturbation. The right panel shows how the Coriolis force makes the fluid rotate in the same direction for both cases, and also that the typical radius of the phenomenon is proportional to 1 / sinϕ.

Fig. 12 Time series of a vertically propagating Alfvén wave. The radial coordinate is expressed as a fraction of one height scale. In the left panel the intensity of the longitudinal magnetic field versus radius is plotted. In the right panel, the longitudinal velocity u is plotted.

Fig. 13 Comparison between the measured Alfvén speed vA from the code (solid line) and the theoretical prediction (dashed line). The time is scaled to the time for the wave to reach the end of the box, since it is difficult to measure the propagation velocity there. It is the same reason why the predicted and measured velocities are not the same at the beginning of the simulation.

6. Conclusions

We present here a numerical magnetohydrodynamic scheme using the hydrostatic approximation, based on the equations developed by Kasahara & Washington (1967) for use in atmospheric circulation modeling. The scheme uses spherical coordinates, and is suited for the modeling of various phenomena in gravitationally stratified fluids. Suitable astrophysical contexts have a strong gravitational stratification and entropy gradient pointing upwards, leading to much greater horizontal length scales than vertical. In addition, timescales of interest should be longer than the vertical acoustic timescale (H_P/c_s), so that hydrostatic equilibrium is established on a timescale shorter than any timescale relevant to the process being studied.

In the magnetic case, a high plasma-β is required for the scheme to be useful. This is because while the hydrostatic scheme filters out vertically propagating sound waves, it does follow vertically propagating Alfvén waves. Consequently, the condition that the time step is limited by horizontally propagating sounds waves rather than vertically propagating Alfvén waves is where Δs and Δr are the horizontal and vertical (radial) grid spacings. Of course, it may still make sense to use this scheme if β is somewhat below this, depending on the parameters involved, but not where β ~ 1.

We have tested the code in various ways. It reproduces phenomena such as the Kelvin-Helmholtz instability, as well as the propagation of a number of different kinds of waves with favorable quantitative comparison to the theory.

This scheme will be useful in stellar interiors, where these stratified conditions apply and where a spherical geometry is useful. This includes conditions close to the stellar surface where approximations such as Boussinesq and anelastic have difficulties, since they are only suitable for modeling phenomena that do not differ much from a constant background state. The hydrostatic approximation can thus be more advantageous since it allows for compressible supersonic flows in the horizontal directions. We are already investigating various (magneto)hydrodynamical processes with this numerical code.

Acknowledgments

We would like to thank Andreas Reisenegger for constructive comments on the manuscript. This work was supported in part by the DFG-CONICYT International Collaboration Grant DFG-06.

References

Appendix A: Interpolations, derivatives, and integrals

We describe in each case two possible scenarios, one in which the variable is cell centered and the result is face centered and the opposite case. Due to our coordinate conventions, the first case will be named “up” and the second “down”. The index convention can be seen in Fig. A.1, capital letters indicate a function defined at the center of the cell and lowercase letters indicate a function centered on the face of the cell.

Fig. A.1 Generic index convention. The coordinate increases in the direction of increasing index.

Appendix A.1: Interpolations and derivatives

Since the coordinates are orthogonal, we discuss only a generic 1D case. Let’s assume we want to know the interpolated value and first derivative of a function f(x) at . We expand using Taylor series around the point we want, using six neighboring points. We get a system of linear equations with the value of the interpolated function and the first five derivatives at the point we want as unknowns, and an error of evaluated at some point of the interval .

We only need the interpolated value and first derivative, so we get two algebraic equations for interpolation, which we call upx (Eq. (A.1)) for and dnx for (Eq. (A.2)), and two for the derivatives dupx (Eq. (A.3)) and ddnx (Eq. (A.4)).

with , , , , and

Appendix A.2: Integrals

The code integrates various quantities in the vertical direction. We will use here a classical approach in numerical integration method. We use the fact that the grid is regular, so we only need to integrate in general in one grid cell. We can expand the function to integrate around the lower integration limit, and then integrate an easy polynomial function. The derivatives and interpolations expressions can be derived from the previous system of equations, so we expand up to sixth order. After simplification and doing some algebra, we realize that the first point is actually not needed, so we end up with operator similar to the previous ones, iupx (Eq. (A.5)) and idnx (Eq. (A.6)). with , and

All Tables

All Figures

	Fig. 1 Positions in each grid cell at which various quantities are defined. The longitudinal, latitudinal and radial components of velocity are denoted by u, v and w.
In the text

	Fig. 2 Contour plot of the pressure in the plane wave test at the center of the radial range. The pressure perturbation moves towards the edges, preserving shape, like a 1D. After one period, the wave is back to its initial position. This behavior is constant in radius.
In the text

	Fig. 3 Fractional difference between the pressure after one period and the initial pressure, taken at a height in the center of the domain. We can see that there is no clear dispersion of the wave packet, as expected from an inviscid fluid.
In the text

	Fig. 4 Contour plot of the pressure in the plane wave test with higher angular range, with the expected theoretical position of a wave traveling with the local sound speed as a black dashed line. The wave appears to move faster away from the equator because the grid longitudinal spacing is smaller by a factor cosϕ.
In the text

	Fig. 5 Contour plot of the pressure in the 2D wave test. We can see how it expands as a circular wave.
In the text

Fig. 6 Shear instability in a square box of 50 × 256 × 256 resolution with a perturbation of n = 5. Top-left: initial longitudinal velocity u. The other three frames are at t = 281, approximately 3.6 sound crossing times: top-right: passive tracer used to follow the locations of the 2 fluids; bottom-left: the longitudinal velocity u and bottom-right: latitudinal velocity v. These horizontal cuts are taken half-way up the computational domain.

In the text

	Fig. 7 Amplitude and growth rate of shear instability modes. Top: amplitude as a function of time, in units of the sound crossing time 0.5τs. Bottom: the growth rate measured after t = 0.6τs (determined by taking the slope of a linear fit of the logarithm of the perturbation amplitudes), are roughly proportional to the wave number n, and roughly agree with the predicted growth rate from Eq. (34).
In the text

	Fig. 8 Shear instability in a square box with the shearing occurring in the radial direction for a model with an initial n = 1–5 perturbation. Panel a) shows the passive scalar at the start of the simulation in an equatorial cut. Panel b) shows the passive tracer after 7.5 sound crossing times, here it is evident that the n=5 mode dominated in the early evolution.
In the text

	Fig. 9 Growth rate of the n = 2–5 modes for the shear instability when the shearing boundary is in the radial direction, for the first 1.8 sound crossing times. The n = 5 mode is the solid blue line, the n = 4 mode is the dashed gold line, the n = 3 mode is the dotted black line, and the n = 2 mode is the dashed-dotted red line.
In the text

	Fig. 10 Pressure normalized to the surrounding pressure versus latitude at the longitude of the perturbation. The southern hemisphere perturbation is a overpressure relatively far from the equator, and the northern perturbation is relatively close to the equator.
In the text

Fig. 11 Streamlines for two snapshots of two rotating over- and under-pressurized features rotating and balanced by the Coriolis force, taken in the middle of the radial domain. The left panel shows how the fluid starts to move away from the center of the southern perturbation, and towards the center of the northern perturbation. The right panel shows how the Coriolis force makes the fluid rotate in the same direction for both cases, and also that the typical radius of the phenomenon is proportional to 1 / sinϕ.

In the text

	Fig. 12 Time series of a vertically propagating Alfvén wave. The radial coordinate is expressed as a fraction of one height scale. In the left panel the intensity of the longitudinal magnetic field versus radius is plotted. In the right panel, the longitudinal velocity u is plotted.
In the text

	Fig. 13 Comparison between the measured Alfvén speed vA from the code (solid line) and the theoretical prediction (dashed line). The time is scaled to the time for the wave to reach the end of the box, since it is difficult to measure the propagation velocity there. It is the same reason why the predicted and measured velocities are not the same at the beginning of the simulation.
In the text

	Fig. A.1 Generic index convention. The coordinate increases in the direction of increasing index.
In the text