The Stagger Code for Accurate and Efficient, Radiation-coupled Magnetohydrodynamic Simulations

We describe the Stagger code for simulations of magnetohydrodynamic ( MHD ) systems. This is a modular code with a variety of physics modules that will let the user run simulations of deep stellar atmospheres, sunspot formation, stellar chromospheres and coronae, proto-stellar disks, star formation from giant molecular clouds, and even galaxy formation. The Stagger code is ef ﬁ ciently and highly parallelizable, enabling such simulations with large ranges of both spatial and temporal scales. We describe the methodology of the code and present the most important of the physics modules, as well as its input and output variables. We show results of a number of standard MHD tests to enable comparison with other, similar codes. In addition, we provide an overview of tests that have been carried out against solar observations, ranging from spectral line shapes, spectral ﬂ ux distribution, limb darkening, intensity and velocity distributions of granulation, to seismic power spectra and the excitation of p -modes. The Stagger code has proven to be a high-ﬁ delity code with a large range of uses.


Introduction
Many astrophysical objects are hydrodynamical in nature, and since hydrodynamics eludes analytical analysis in all but the simplest cases, simulations are indispensable for the analysis and interpretation of observations of these.We present here a code for such simulations, the Stagger code, named after the staggering between dynamical and thermodynamical variables in the simulation domain.This code is modular, highly parallelizable, has accurate magnetohydrodynamic (MHD) solvers, and is adaptable to many kinds of astrophysical systems.
The Stagger code has also taken part in code comparison and verification efforts, most notably by Beeck et al. (2012) and Kritsuk et al. (2011).
The Stagger code has previously been described briefly by Nordlund & Galsgaard (1995), with the present document being a comprehensive update, aiming for a more complete presentation.The structure and organization of the code is explained in Section 2, the governing equations are laid out in Section 3, spatial operations are presented in Section 4, how time is advanced is described in Section 5, the options for artificial viscosity are presented in Section 6, the physical heating and cooling mechanism are described in Section 7, the available microphysics is outlined in Section 8, external and internal forcing options are listed in Section 9, boundary conditions are discussed in Section 10, file handling is explained in Section 11, a large range of numerical and observational tests are analyzed in Section 12, and we provide some perspectives in Section 13.

Structure and Philosophy of the Code
The Stagger code solves the MHD equations.Its basic philosophy is efficient parallelization and modularity.This involves using the Message Passing Interface (MPI) and segregating different functionality-and environment-dependent aspects.The goal was to enable its easy use to simulate a wide variety of scenarios on computers ranging from laptops to massively parallel machines.
All operating-system-dependent aspects are segregated into files containing any system-dependent Fortran functions, appropriate compiler commands and flags, and system commands.MPI routines are called only from the differentiation, extrapolation, and I/O routines, which are contained in separate subroutines and directories.The MPI routines themselves are in a separate directory.The code to calculate the time derivatives of the physical variables is in two subroutines, one for the fluid and the other for the magnetic field.Time-stepping procedures, equation of state (EOS) applications, heating and cooling, and external forcing are all done in separate subroutines.Each experiment has a subdirectory with its own specialized routines, including any specialized boundary conditions.
The Stagger code routines are organized in a tree.At the top are the basic routines common to all calculations-those that calculate the time derivatives, do the time stepping, do the I/O and the MPI operations.At the next level are separate directories for the spatial derivatives and interpolation, which call MPI operations; the EOSs; heating/cooling; external forcing, e.g., gravity; initial set up; utilities such as fast Fourier transforms; and finally a tree of set-up instructions and specific routines for different astrophysical experiments, e.g., magnetoconvection, star formation, and galactic structure.
This allows the Stagger code to be applied to many different situations with different physics, different geometry or computational mesh, and for it to be run on a wide variety of computational platforms, with a need to change only a few routines.It can be run on most computers, from PCs to massively parallel supercomputers.The modularity also makes it easier to maintain and debug the code, as problems can be localized and addressed quickly, considering the ∼280k lines of code involved.The Stagger code is maintained in a git repository git@bitbucket.org:aanordlund/stagger.git.The version used for testing is tagged 2024-ApJ.

Equations
The Stagger code calculates the nonsplit time derivatives from the conservation equations for mass, momentum, internal energy, and the induction equation for the magnetic field, together with Ohm's law for the electric field.
Mass conservation controls the topology of stratified convection, where ρ is the density and u the velocity.Momentum conservation controls the plasma motions: Here P is the pressure, g is the gravitational acceleration, B is the magnetic field, J = ∇ × B/μ is the current, μ is the permeability (magnetic constant), 2ρΩ × u is the Coriolis force, and is the viscous stress tensor, is the strain tensor, ν ij the coefficient of viscosity, and c 1 is an input parameter.Typically, we take c 1 = 0.If no viscous pressure resistance to overall contraction is desired, then c 1 = −1/3.In most of the computational domain, we use sixth-order derivatives in the viscous stress tensor.However, where the magnitude of the second-order differences is large, we switch to second-order derivatives.
At large depths where fluctuations are small there can be issues with round-off error in the derivatives of the pressure.To address this, Stagger has a switch to do all the derivative and centering operations on ( ) P P ln aver , where P aver is the horizontally averaged gas pressure at the given vertical position.For instance, the vertical derivative of the pressure becomes The horizontal derivatives are similar except that the horizontal derivative of P aver vanishes.Internal energy is changed by transport, by PdV work, by Joule heating, by viscous dissipation, and by heating and cooling.It is the fluid version of the second law of thermodynamics and (together with the density) determines the plasma temperature, pressure, and entropy: where e is the internal energy per unit mass, and Q is the heating or cooling.This can be Newtonian cooling or by thermal conduction, optically thin radiation, or by radiative transfer.
The viscous dissipation is Convection influences the magnetic field via the curl E term in the induction equation: the electric field is given by Ohm's law.In a one-fluid MHD system, it is and η is the resistivity, n e is the electron number density, P e is the electron pressure, and e is the electron charge.The last two (Hall and pressure) terms are usually neglected, but the Hall term may be important in the weakly ionized photosphere.B is periodically cleaned to keep ∇ • B = 0 (Brackbill & Barnes 1980;Tóth 2000).Since B = ∇ × A + ∇f (where A is the vector potential and f is the scalar potential for the magnetic field), This is solved for f in Fourier space, and the corrected B is In addition, a term can be added to the magnetic field time derivative,  , (where C divb is an input constant <0.11 for stability, dx, dy, dz are the local grid spacings in their respective directions, and dt the time step) to diffuse away the divergence of B and ensure it is zero to machine accuracy at all times.
The time derivatives of density, momenta, and internal energy are calculated in one routine and the time derivative of the magnetic field in another to enable both fluid and MHD calculations to be easily performed.

Centering
The physical variables are staggered.Density, energy, and pressure are located at cell centers.Momenta, velocities, and magnetic fields are centered on cell faces at the lower index side.Electric fields are centered on the cell edges (Figure 1).

Derivatives
Derivatives, except possibly in diffusion terms, are evaluated by sixth-order, centered differences in index space with the Jacobian of the transformation back to physical space: The weighting coefficients are as follows: Similar formulas are used for each dimension and for downward derivatives.The Jacobian for equally spaced directions is simply 1/δx i .For stretched grid directions, the Jacobian is calculated by first making a spline fit to the grid spacing and taking the Jacobian as the inverse of the fitted grid interval at the center of the derivative.All the derivative routines assume periodicity and take variable values falling outside the boundary from the corresponding locations in the opposite boundary.MPI routines are used get these values.In nonperiodic directions, boundary conditions load values into ghost zones at each boundary.

Interpolations
Interpolations to properly center variables are calculated by a fifth-order, centered template in index space:

Coordinate System
In stratified calculations (e.g., stellar convection) the stratified direction is taken as the y-direction.The x-and zdirections are periodic.In stellar calculations the indices increase downward into the star.Hence, the lower boundary in index space is in the photosphere and the upper boundary deep in the convection zone.Together with the horizontal xand z-directions, this forms a right-handed coordinate system.

Time Advance
Time advance is by a low-memory-usage, third-order Runge-Kutta scheme (Kennedy et al. 1999).The time derivatives of the variables are evaluated in a nonsplit manner, three times for each full time step, using Equations (1)-( 8).The time advance scheme is as follows.Suppose the time derivatives are given by a function F, i.e., In a cycle of three iterations, i = 1,3, updated derivatives are calculated and added to a fraction, α i , of the previous derivatives and stored in the same location: are actually the same location in memory.First, the existing time derivative is multiplied by α i and restored in the same location.Then, each term in the new time derivative is added to ( ) f t i ¶ ¶ as they are calculated.In the first iteration α 1 = 0, so for each time step all the derivatives are first initialized to zero.After the update of the time derivatives is completed, the variables f are updated: and the next iteration begun until after three cycles the variables for the next time step are obtained.
Only one set of variables and one set of their time derivatives needs to be stored.

Courant Condition
The time step is set by a combination of the usual Courant condition with the wave speed plus fluid velocity (in the case of magnetic field the fluid velocity perpendicular to the field), together with a limit depending on the rate of change of the density and energy and a limit depending on the rate of dissipation (Equation ( 19)): Here min ijk and max ijk are the minima and maxima over the entire computational volume (excluding ghost zones).C dt , C dtr , and C dtd are input parameters (see Appendix A.2). C s is the sound speed and C A is the Alfvén speed.ν and η are the viscosity and resistivity.

Viscosity
Higher-order schemes (higher than linear) are inherently unstable, and need a mechanism to suppress the growth of unphysical features.One such mechanism is an artificial viscosity, to dampen unphysical oscillations between neighboring grid points, ringing at sharp changes, and to stabilize sound and Alfvèn waves as well as shocks.
There are two philosophically different approaches to artificial viscosity.One is for the viscosity to describe the effects at the grid scale, of the actual, molecular viscosity acting on the dissipation scale of the plasma (see Section 6.1.2).Another is the entirely utilitarian approach of stabilizing the simulation with a minimal amount of viscosity, of which Section 6.1.1 gives an example.
There are input variable switches to choose different types of viscosity (see Appendix A.5).The simplest is a constant viscosity const ij n = .This is useful for testing purposes or for making direct numerical simulations, but is otherwise an inefficient use of a given number of grid points in a simulation.In general, we use a Richtmyer & Morton (1967) type of viscosity proportional to the wave and fluid speeds.The viscosity model by Richtmyer & Morton (1967) is where s is the sound speed and a is the Alfvèn speed, and where c v , c sh , and c w are input parameters (nu1, nu2, and nu3).
The viscosity is scaled by the grid spacing, Δx i , which accommodates nonisotropic grids, and also makes the choice of diffusion constants, e.g., nu1, nu2, nu3, independent of resolution to first order.The influence of the local viscosity is spread by taking its maximum over three or five adjacent cells and then smoothing it by taking the average over three adjacent cells in each direction, denoted by : For the off-diagonal viscosity, the averaging is carried out on ln n and then the exponential is taken, in order to avoid edge effects of large gradients.

Smagorinsky-type Viscosity
Another possible form for the viscosity is the Smagorinsky (1963) type.We take This viscosity is sometimes reduced with depth in stratified atmospheres by a factor [ , where ρ 0 is an input reference density.In this case, also, the influence of the viscosity is spread by taking the maximum over three or five adjacent cells and then smoothed by averaging over three adjacent cells: Application of the Smagorinsky approach is often justified by it describing the viscosity in the case of a self-similar turbulent cascade of isotropic turbulence.For a numerical simulation to exhibit isotropy at the grid scale, it would need a large inertial range, though making this argument spurious in many cases.This does not, however, affect its utility in stabilizing a numerical scheme with little effect in smooth regions.
We often include both the Richtmyer & Morton and the Smagorinsky viscosities in the same simulation in varying proportions.

Quenching
We reduce the viscous stresses in smooth regions (and enhance it for two zone wiggles) by the ratio of the second differences to the absolute value of each viscous stress tensor component, T: q T T T q T q q q q q min max , , , , 24 where q max and q lim are input parameters (see Appendix A.8).

Diffusion in the Energy Equation
The momentum (Equation ( 2)) is stabilized by dissipation (Equation (3)), a diffusive process.This dissipation acts as a heat source in the energy (Equation ( 5)).To keep the energy Equation (5) stable requires the addition of some diffusion, x e e x ln , 2 5 where c E is an input parameter (nuE in Appendix A.5).When there is a mean gradient in the energy, we want to only diffuse fluctuations about the mean.In that case, the energy diffusion is x e e e x ln ln , 26 where e ln horiz á ñ is the horizontal average of the internal energy per unit mass.

Newtonian Cooling
Newtonian heating/cooling makes the energy tend toward a fixed value e cool on a timescale t cool : Newton cool cool

= --
This can be multiplied by a function of density to smoothly reduce the cooling at either low or high densities, ρ cool , e.g., , where the cooling time t cool defaults to 0.1 s.

Thermal Conduction
Heating/cooling by thermal conduction is given by The conductive flux, F cond , in the presence of a magnetic field, B, is where B  is the unit vector in the direction of the magnetic field.Spitzer conduction (Schrijver & Zwaan 2000) for a dilute plasma is implemented as where k is Boltzmann's constant, m e the electron mass, e the charge of an electron, and the temperature is approximated as 32  gand, for this particular case, the adiabatic exponent is approximated as a constant, γ = 5/3.

Optically Thin Cooling
If radiation can easily escape from the entire volume, then its only effects is to cool the plasma, in a way that depends exclusively on local conditions: (where Λ is a function of the local density and energy).This is appropriate in, for example, the solar corona and the interstellar medium.For the corona, the formula used is ´----- (Kahn 1976), where the number density of all atoms and ions for a representative composition is approximated with ∑ i n i = ρ/(2 × 10 −24 g).The exponential quenching factor suppresses the cooling for T < 10 4 K, to better approximate the basis for Kahnʼs (1976) original expression without it.The interstellar medium is another milieu in which optically thin cooling by various metal atoms and ions is appropriate.In this case, the cooling function is more complicated, with different physical processes being important in different temperature regions; see Dalgarno & McCray (1972).

Radiative Transfer
When neither the optically thin nor the optically thick diffusion limit is appropriate, the radiative heating/cooling is where I λ is the monochromatic radiation intensity in a beam dΩ in a given direction, J λ is the mean intensity averaged over all directions, and B λ is the Planck function.The total monochromatic opacity, ρκ λ = a λ + σ λ , is composed of true, thermalizing absorption, a λ , and isotropic and coherent (no change of photon wavelength) scattering, σ λ .Both absorption and scattering coefficients are per volume, but the total extinction, κ, is per mass, hence the factor of density.
The radiative transfer equation describing I λ is where ℓ is the path length along the ray.One can divide by the opacity ρκ λ to get the transfer equation in terms of optical depth, which is the form used to solve it, where dτ λ = ρκ λ dℓ is the optical depth (the distance in units of the photon mean-free path).μ and f specify the direction of the ray along which the transfer occurs.Here, is the probability a photon will be destroyed (i.e., thermalized).So far, the Stagger code assumes local thermodynamic equilibrium (LTE), which means ò λ = 1, but the option of ò λ ≠ 1 is left open in the radiative transfer tables of Section 8.2.A non-LTE treatment was devised by Hayek et al. (2010) for the Stagger-code-derived Bifrost code (Gudiksen et al. 2011).
The radiative transfer employed by the Stagger code has been described in detail (e.g, Nordlund 1982;Nordlund & Dravins 1990) for earlier simulations.The forward solution is performed on long characteristics, using a Feautrier technique (Mihalas 1978), modified to give more accurate solutions in the optically deep layers (Nordlund 1982).The monochromatic solution is prohibitively expensive for MHD calculations involving thousands of time steps.The Stagger code speeds up the solution enormously by reducing the many wavelengths used in the solution by grouping them into bins of similar opacity (i.e., opacity binning or a multigroup method) with only a few bins and by approximating the 3D angular distribution with a vertical ray and only a few azimuthal rays.
The angular integration in the azimuthal direction, j, is carried out with a rotating set of equidistant angles in order to minimize directional biases.The integration over inclination angle cos m q = (μ = 1 is the vertical ray) is performed with Radau quadrature (Radau 1880;Abramowitz & Stegun 1972), which includes the vertical ray and N μ − 1 slanted rays.It is implemented for N μ = 1-10.A common choice is N μ = 3, N j = 4.A quadrature scheme by Carlson (1963) and Bruls et al. (1999) is also available for N μ = 2-4.
Radiative transfer is solved in detail for layers where ( ) min 500 horiz max , and ignored in layers deeper than that.In the deeper layers the diffusion approximation is valid, but radiative heating and cooling is insignificant there for cooler stars.For stars hotter than T eff ; 7000 K the radiative heating/cooling in the diffusion approximation can be included as the divergence of the radiative flux: where here σ is the Stefan-Boltzmann constant and κ Ross is the Rosseland mean opacity: Being a harmonic mean, the Rosseland opacity is heavily weighted toward the lowest opacities, typically located in the optical part of the spectrum.The Planck mean, on the other hand, is a straight average:

. General Considerations of Opacity Binning
Opacity binning is a method for dramatically reducing the computational cost of realistic radiative transfer from that of the monochromatic case with of order 10 5 wavelength points to a binned case with a mere dozen bins, N bin .Each bin, labeled i, has a weight, w i , of dimension wavelength.
The bin-wise Planck functions, B i , and bin-wise weights, w i , are computed as where {λ} i denotes the subset of wavelength points that belong to bin i, and w λ is the weight of wavelength λ.The total Planck function is a simple sum of the bin-wise contributions: as is also the case for the total specific intensity I, angularaveraged intensity J, flux F, and heating Q.These latter four quantities are all the result of radiative transfer calculations on the bin-wise B i and bin-wise τ i scales, and the τ i in turn results from the integrations of the bin-wise opacities κ i .Different binning schemes are defined by how bin membership is decided for each wavelength, and how bin-wise opacities and source functions are computed.These are based on a monochromatic reference calculation (see Section 7.6), and the binning schemes aim to reproduce the cooling and heating of this reference calculation.The Stagger code currently offers the choice of two schemes, as detailed in Sections 7.7.1 and 7.7.2.

Monochromatic Reference Calculation
The wavelengths are sorted into bins based on a monochromatic, one-dimensional (1D), reference radiative transfer calculation, performed on the averaged atmosphere structure of a simulation.The average structure is a horizontal mean of 〈ρ〉 and 〈T〉, averaged in time to reduce p-mode effects.The temporal averages are performed on a horizontally averaged column mass scale and interpolated back to the geometric scale.This is a pseudo-Lagrangian temporal average that follows the mean vertical motion, effectively filtering out firstorder effects of the p-modes.This reference model is stored for later use, and only updated if the atmosphere of the simulation has drifted, as will be the case during relaxation of a new simulation.During production runs, however, the reference model should not be changed.This makes the tables of binned opacities specific for a given simulation, and it makes little sense to use it for another set of atmospheric parameters.
The 1D radiative transfer calculation on the reference model results in two-dimensional (2D) arrays (as a function of wavelength and depth) of monochromatic absorption coefficients a 1D l , scattering coefficients 1D s l , optical depths 1D t l , and zeroth and second angular moments J 1D l and K 1D l , respectively, of the specific intensities I 1D l .This radiative transfer calculation is carried out for N μ = 4 Radau (1880)-distributed μ angles, and since it is in 1D no explicit f integration is needed.

Sorting Wavelengths into Bins
A particular wavelength is ascribed a bin according to its monochromatic opacity strength.This opacity strength, x, is based on results of the 1D reference calculation of Section 7.6, and evaluated as either the monochromatic extinction in the monochromatic photosphere at wavelength λ.
The Stagger code gives the option between two methods of opacity binning: binning with scaled opacities as described by Nordlund (1982) and in Section 7.7.1, using Equation (45), and binning with individual opacities, which is based largely on the work by Skartlien (2000) but with additional improvements as described in Section 7.7.2 and using Equation (44) together with bin divisions in wavelength.7.7.1.Opacity Binning: Scaled Opacities Nordlund (1982) introduced the method of opacity binning as a first attempt at a realistic treatment of line blanketing in 3D convection simulations.The complete spectrum is binned according to the monochromatic opacity strength, x log 10 l , of Equation (45), based on an optical depth of the Rosseland mean across all wavelengths, but excluding line opacity.The bins are assumed to be equidistant in x log : 10 ( ) x i 10 , with 1 and 0, 1, 2, 3.
Furthermore, the bin-wise opacities are approximated as simply proportional to the i = 0 (continuum bin) opacity: This was done to optimize the convection code for speed, memory, and disk space, in order to make the original convection simulations (Dravins et al. 1981;Nordlund 1982;Nordlund & Dravins 1990) feasible at all.The subsequent increase in computational power has made it possible to abandon these simplifications, as described in Section 7.7.2,below.
The continuum bin opacity, κ 0 , is computed as a smooth transition between a Rosseland mean of the continuum, 0  k , below and an intensity mean, 〈κ〉, above the photosphere: where 0  t is the optical depth integrated from 0  k of Equation (49).This ensures correct asymptotic behavior in the free-streaming regime in the optically thin case, as well as in the diffusion regime in the optically thick case.The bridging constant, c τ , has a limited effect on the results and defaults to a value of c τ = 2, which has been tested against 1D reference models for a range of stellar atmosphere parameters.
The bin-wise Rosseland opacity is This definition is not confined to the domain of the 1D reference calculation, since we can evaluate B i for any temperature and therefore for every point in a 2D table.
For the scaled opacity binning method, only the continuum mean 0  k is used, resulting in the effective bin-wise opacity, Equation (47) via Equation (48).
The intensity mean, 〈κ〉, is an average over all wavelengths: where J 1D l is the angular average of the monochromatic specific intensity, I λ .The exponential factor ensures that the contribution of 〈κ〉 to κ 0 is suppressed in the photosphere and below.J 1D l is a result of the 1D reference calculation of Section 7.6, and hence necessitates an extrapolation away from that 1D stratification to span the 2D table, as detailed in Section 7.7.3.It is not, however, 〈κ〉 1D itself that we extrapolate, but rather a correction factor: from which the total, bin-wise opacity can be computed as and where x corr is the extrapolation of x corr 1D away from the 1D reference model.
The harmonic mean of the bin-wise, scaled 0  k Rosseland opacity does not in general converge to the total κ Ross in the limit of many bins, but rather overestimates κ Ross from the photosphere and below.This overestimation can range from factors 2-5 in the photosphere, depending on stellar atmospheric parameters.

Opacity Binning: Individual Opacities
This method generalizes the bin assignments to also include wavelength.This was done in recognition of various prominent bound-free absorption edges in the ultraviolet, which can give rather different behavior for opacities between the blue and red sides.The molecular bands in the infrared also exhibit different behavior than atomic lines in the optical.An example of this generalized binning scheme is shown in Figure 2 for a solar simulation.The case of Equation (46) can be trivially reproduced in this more general version.The binning can be optimized, i.e., moving the borders between the bins in order for the radiative cooling of the binned solution, to reproduce that of the full monochromatic solution.In Figure 3, we show an optimized solution, corresponding to the binning shown in Figure 2.
The cooling peak, seen in Figure 3, shapes the photospheric transition.It is crucial to reproduce this feature accurately in order for the simulation to equilibrate to the correct T eff and in order to produce realistic synthetic spectra from the simulation.2. The differences between the monochromatic and the binned solutions are shown, exaggerated by a factor of 10, with the red and blue dashed lines on the right-hand scale.Left: the four equidistant bins of Section 7.7.1, also shown in the left-hand panel of Figure 2. Scaled opacity calculation shown in blue, and individual bin-wise opacities, as in Section 7.7.2, in red.With scaled opacities the maximal difference with respect to the monochromatic calculation, relative to the magnitude of the cooling peak (shown with vertical dotted line), is 21.75%, the overall rms deviation is 6.58%, and the rms deviation just in the cooling peak is 10.32%.For individual opacities the deviations are 2.81%, 0.63%, and 1.13%.Right: a 12-bin implementation of the optimized binning of Section 7.7.2, also shown in the righthand panel of Figure 2, giving deviations of just 1.22%, 0.26%, and 0.45%, respectively.
For the intensity mean this scheme uses the bin-wise version of Equation (50): This opacity does not include scattering, unless the table is being made for non-LTE source functions.For that case, scattering is simply added to the absorption above.
As for 〈κ〉, the bin-wise opacity is also confined to the domain of the 1D reference calculation, and an extrapolation is needed in order to span the extent of the table.
Analogous to the case of scaled opacity binning, the quantity we extrapolate is the correction factor, x i corr, 1D , that turns the binwise Rosseland mean, i  k , into the effective bin-wise opacity, κ i : The extrapolation of this to span the 2D plane of the table is described in Section 7.7.3.The final bin-wise opacity is In contrast to the binning with scaled opacities, this version of κ i converges properly to the Rosseland mean in the limit of N bin → ∞ .This ensures the photosphere emerges in the correct location, whereas with scaled opacities it is too high up in the atmosphere, i.e., at lower pressure/density.This version with individual opacities was employed by, for example, Asplund et al. (2009), Collet et al. (2011a), and Magic et al. (2013), and is a major cause for the differences between the solar abundances of Asplund et al. (2005) and Asplund et al. (2009).Collet et al. (2018) performed a thorough analysis of different layouts and numbers of bins, concluding that even 48 bins does not result in radiative heating that is significantly closer to the monochromatic calculation, compared to the standard 12 bins.

Extrapolation of the x corr 1D Factor
We need to extrapolate the factor x corr 1D from its 1D stratification and out to the 2D (ρ, e) plane of a table, so that it can affect the whole simulation in a plausible way.In Figure 4, we show the distribution of ρ and T for 17 horizontal planes in a convection simulation.It is clear that the distribution of points is far from uniform, and we can fit them fairly well with which connects a point (T, ρ) in the table to a point j in the reference stratification, via the fitted black lines of Figure 4 and the coefficients a 1 -a 5 .

Vertical Scale for Radiative Transfer
The radiative transfer calculation can often benefit greatly from a vertical scale that is explicitly optimized for that purpose.There is therefore an option for stellar atmosphere calculations to compute such a vertical scale for the transfer calculations that better resolves the radiation field (Collet et al. 2018).The radiative heating/cooling is evaluated on this dynamically optimized radiative grid, and then interpolated back to the hydrodynamic grid of the simulation for addition to the energy conservation (Equation ( 5)).

Equation of State and Opacity
Simple opacities, such as H − , can be calculated analytically as a function of ρ and T as needed.Complete stellar atmosphere opacities, however, are too complex for that and instead auxiliary codes compute tables which are used by the Stagger code for looking up opacities and EOS variables.
The Stagger code also has some simplified options for EOSs with only a few elements that are computationally cheap to solve.
A couple of more realistic sets of atomic physics, relevant for stellar atmospheres in particular, are elaborated on below.

Original Tables for Stellar Atmospheres
The original version of atomic physics tables for stellar atmosphere simulations consists of a table in ln r and θ = 5039.78/Tthat gets inverted into tables in ln r and ρe as needed by the Stagger code.Both tables contain gas pressure p ln , continuum bin opacity ln 0 k , the N bin source functions B ln i and the complement to the independent temperature variable ρe and T, respectively.Each variable in the final table is also accompanied by numerical d ln r and d(ρe) derivatives, to enable bicubic interpolation in the tables by the Stagger code.Only the final table is used by the Stagger code, but the initial table can be useful in postprocessing.
These tables are specific to a particular simulation, and have irregular shapes that just covers the extent of the simulation, with a small buffer.If the simulation strays outside this table it is automatically expanded where needed.The θ-resolution of each isochore in the table is determined from a requirement that θ-interpolations in any of the table variables do not deviate more than ò from a direct calculation.This parameter is typically set to ò = 10 −5 .
This version assumes the scaled opacity binning of Section 7.7.1, with LTE source functions, S i = B i .Both the opacities and thermodynamics are based on the atomic physics of the original MARCS stellar atmosphere package (Gustafsson 1973;Gustafsson et al. 1975).
EOS.This EOS is based on an interpretation of the Debye & Hückel (1923) theory of electrolytes, which lowers the ionization potential of all species according to the formulation by Griem (1964).The partition functions of atoms and ions include an asymptotic part that includes some accounting for the abovementioned potential lowering (de Jager & Neven 1960).This alters the populations of the various species, but there are no direct effects on the internal energy and the pressure, making this EOS thermodynamically inconsistent (e.g., values of the adiabatic gradient will depend on which valid expression is used).
These tables are computed for a 16-element mixture of H, He, C, N, O, Ne, Na Mg, Al, Si, S, K, Ca, Cr, Fe, and Ni.Molecules of the most abundant elements, H 2 , H 2 + , H 2 O, OH, CH, CO, CN, C 2 , O 2 , N 2 , NH, and NO, are included by means of the equilibrium constants by Tsuji (1973).
Opacities.The bound-free (bf) opacities are listed in input tables as a function of wavelength, and also as a function of T if the source has excited states.The included sources are H 2 + bf and free-free (ff), H 2 -ff, C I, Mg I, Al I and Si I bf, H − bf and ff, as well as the free-free and bound-free opacities of the first 15 states of hydrogen.Scattering by free electrons and Rayleigh scattering by H and H 2 is included via analytical expressions.Spectral lines are included in the statistical approach of opacity distribution functions (ODFs), with 42 distribution functions spanning 3000-7200 Å, each 100 Å wide and sampled at four points.These ODFs are tabulated for nine temperatures T = 3-9 kK and 15 electron pressures p log e = 3 --4.They involve about 50,000 atomic lines on more than that of molecular lines of MgH, CH, OH, NH, and CN.For λ = 7200-125000 Å, a separate set of 23 ODFs are used for lines of CN and CO molecules, only, provided for six temperatures between 2 and 7 kK.
The opacity in these tables is the continuum bin opacity, κ 0 , of Equation (52).

New Tables for Stellar Atmospheres
This version consists of a rectangular table in ln r and T ln containing quantities needed for radiative transfer (source functions and opacities for each bin or wavelength), and a rectangular table in ln r and ( ) e ln r of the EOS, as well as some opacities for postprocessing.
Tables have been computed for a variety of scaled solar abundances of the 17 elements H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, K, Ca, Cr, Fe, and Ni by Trampedach et al. (2014).The first simulation to be computed with these tables is the solar atmosphere that the Asplund et al.EOS.The EOS table contains Ntbvar 6 rectangular arrays of gas pressure P ln gas , Rosseland mean opacity ln Ross k , temperature T ln , electron number density n e , Planck mean opacity ln Planck k , and the monochromatic opacity at 5000 Å, ln 5000 k .Other variables can be added at the end.These variables are functions of a regular grid of density, ln r, and internal energy per mass, e ln .The extent and resolution of the grid is given in the header of the table, together with the number of elements and their abundances used in the calculation, the option flags that were set for the calculation, and the revisions of the codes used in the calculation, as well as date and timestamps, uniquely identifying the table.
The thermodynamics is supplied by the Mihalas-Hummer-Däppen (MHD) EOS (Hummer & Mihalas 1988;Mihalas et al. 1988), which is a model for the Helmholtz free energy, from which all thermodynamic quantities can be evaluated consistently via the Maxwell relations (e.g., Däppen et al. 1988).Nonideal effects (a.k.a.pressure ionization) is physically and smoothly accomplished through the so-called occupation probabilities assigned for each excited state, and based on a formulation of perturbations of bound levels from the fluctuating electric field of passing ions.This gives both a physical and differentiable truncation of the otherwise diverging partition functions.All ionization stages and all excited levels are explicitly included, but only two molecules, H 2 and The points are a random subset of the simulation cube&#39;s upper 17 layers, and the lines are the fits described in the text, with coefficients a 1 -a 5 .This example is for the simulation t47g22m+00 at t = 10 6 s.
H 2 + , are treated in the MHD EOS.The independent variables are T log 10 and log 10 r, which need to be inverted to ln r and e ln for the Stagger simulations.This is done by exploiting the analytical partial derivatives of the MHD EOS, which allows for an accurate and physical third-order interpolation to a new regular e ln grid, and with analytical and consistent e ln derivatives.
Convective stability, particularly in deep layers where the velocities and entropy fluctuations are very small, are sensitive to errors in the EOS, so a bicubic interpolation in the tables is used.Each  , and a resolution of N e = 300 and N ρ = 57.The latter is dictated by the original MHD EOS calculation, which typically uses Δρ = 0.25.
Opacities.The tables are rectangular and regular in the independent variables ln r and T ln .They contain the relative opacity, x i = κ i /κ Ross , of bin i, the Planck part of the source function, ( )  B ln B i i , , with B i from Equation (42), and the photon destruction probability, ( )  ln J i , (Equation ( 38)), of each of the N bin bins.The tabulation of ò B,i allows for running simulations including scattering, although this is not yet implemented in the Stagger code.The LTE case is recovered when all ò B,i = 1 (σ i = 0).The three last entries are occupied by the opacities ρκ Ross (Equation ( 40 Since the evaluation of radiative transfer in an atmosphere simulation takes up a significant part of the computational cost, and interpolation in the opacity table is a significant part of that, it is most efficient to increase the resolution (compared to that of the EOS table) so only linear interpolations in the table are needed (and no derivatives are stored in these tables).The simulations are not critically sensitive to small interpolation errors in opacity.
The extent and resolution of a table is stored in its header, as are the elemental abundances, the option flags for the calculation, the revisions of the codes used in the calculation, as well as date and timestamps, uniquely identifying the table.
Tables are calculated prior to the simulation runs.Line opacities are currently supplied by the MARCS stellar atmosphere modeling group (Gustafsson et al. 2008) in the form of opacity sampling tables on a 108,405-point wavelength grid covering from 910 Å to 20 μ equidistantly in log l.The 58 included sources of continuum opacities are added in, as detailed in Hayek et al. (2010), and includes molecular sources and collision-induced absorption.The photoionization cross sections of metals have been sourced from the NORAD database by Nahar (2011), based on close-coupling (CC) calculations, abandoning the usual LS-coupling approximation.
Since we want to include molecules in the opacities, and only H 2 and H 2 + are included in the MHD EOS, we revert to the original EOS from Section 8.1 to compute the population of ionization and dissociation stages, based on the electron pressure from the MHD EOS.This introduces a slight inconsistency in the atomic physics.

Free EOS and Blue Opacities
Recently, Zhou et al. (2023) implemented the Free EOS by Irwin (2012), which is constructed to emulate various commonly used EOSs, including the MHD EOS described in Section 8.2.The Free EOS is simpler and faster than the MHD EOS, and the code is publicly available, enabling calculation of tables for custom abundances.It also includes five more elements than the custom-computed MHD EOS tables of Section 8.2, but still only the two molecules H 2 and H 2 + , and no H − anion.This necessitates a transition to an atmospheric EOS for T < 10,000 K that does include a comprehensive set of molecules, as well as H − .Opacities and atmospheric EOS are provided by Blue (Amarsi et al. 2016a(Amarsi et al. , 2018)), which is a major update to the atomic physics package employed by the MARCS 1D stellar atmosphere code (Gustafsson et al. 2008).Blue includes opacities from 21 diatomic molecules and H 2 O, and mostly the same continuum opacities as described in Section 8.2.Atomic and ionic lines are an improved implementation of the same line data also used for both MARCS and the Stagger code physics of Section 8.2.
This package is a standalone code that produces tables in the Stagger table format, and is not part of the Stagger repository.Please address inquiries about this package to its authors.

Forcing
Forcing physics is specified by a separate subroutine.Many routines already exist in the directory FORCING (Table 1).
Additional routines can be written for other specific situations.

Boundary Conditions
The Stagger code has been used for many different types of simulations, each with its own set of boundary conditions.Variable values extending beyond the boundary are needed by the derivative and interpolation routines.The x-and zboundaries are always periodic.The y-boundary can be periodic, closed, open in various ways, or forced.Boundary conditions, when needed, are imposed in the y-direction by specifying variable values in ghost zones at each boundary.The number of ghost zones is determined by the number of derivative and interpolation operations on a variable that are needed for each time step, as each operation destroys data in the outermost valid ghost zone.If the domain is stratified, the stratification is taken to be in the y-direction.The MHD equations with the sixth-order derivatives and fifth-order interpolations require five ghost zones.In this case, the actual boundary point is five zones in from the edge of the numerical box.The so-called lower boundary is located at small indices in the y-direction, i y = lb, with ghost zones i y = [1; lb], and is normally at the top of the simulation domain, as determined by the direction of gravity, if applied.

Periodic
If no boundary conditions are imposed, then all three directions are periodic and there are no ghost zones.The needed variable values for derivatives and interpolations that extend beyond a boundary are taken from the corresponding locations on the opposite boundary.This is the case, for example, for simulations of turbulence, star formation, and disks.

Closed
To impose a closed boundary, the density or its log is extrapolated with a fifth-order extrapolation.The velocity components parallel to the boundary are symmetric about the boundary point and the normal component of the velocity is antisymmetric about the boundary.The energy per unit mass is extrapolated and the energy per unit volume is obtained by multiplying by the extrapolated density.The parallel components of the magnetic field are antisymmetric about the boundary point.

Rigid
For a rigid boundary, the velocity in the ghost zones is taken to be the velocity at the boundary point, and the time derivatives of all variables in the ghost zones are set to zero.

Stellar Atmosphere
Stellar atmospheres are highly stratified.The direction of gravity is taken to be downward in the y-direction of increasing index.In this case, boundary conditions, via ghost zones, must be imposed in the y-direction.Time derivatives of all the variables are set to zero in the ghost zones.

Interior Boundary inside the Convection Zone
Density and energy are needed to relate the momenta and velocities and to determine the pressure in the ghost zones.Their ghost-zone values are determined by linear extrapolation of the log density and the energy per unit mass.The momenta are calculated to be symmetric about the boundary point and velocities are calculated from the momenta and density.
In the inflows, the density and energy are specified (e bot and ρ bot ) in order to fix the pressure and entropy: Outflows should be as unconstrained as possible.However, horizontal pressure fluctuations need to be minimized for stability, without affecting their entropy.This is done as follows.Entropy where E is the internal energy per unit mass and e the internal energy per unit volume.Thus, perturbations at constant entropy must have   the enthalpy per unit mass.We also want to cancel the horizontal fluctuations in pressure in the outflows: ( )  that needs to be added to the existing derivative at the boundary can be obtained using Equation (59): These changes preserve the boundary entropy of the fluid.
The magnetic field time derivative is the curl of the electric field, which preserves ∇ • B = 0.However, we often want to impose a boundary condition on the magnetic field.We therefore store the original magnetic field values in the ghost zones and apply our boundary conditions to the magnetic field.First, extrapolate all three magnetic field components into the ghost zone.Then, set the horizontal magnetic field components to their prescribed boundary values B bdry (Bx0, By0, and Bz0; see Appendix A.6):   et al. 2018).The good agreement shows that the thermodynamic structure of the simulation near the surface is realistic.
where f is a factor that transitions smoothly from zero in the outflows and 1 in the inflows: The electric field in the ghost zones is then given by Ohm's law from the velocity and magnetic field ghost-zone values.
Alternatively, the electric field can be calculated by solving The original ghost-zone values of the magnetic field are restored and ∂B/∂t is calculated from ∇ × E.

Exterior Boundary in the Atmosphere
The log density and internal energy per unit mass are extrapolated into the ghost zones.This extrapolation may be unstable, so the values of the ghost-zone density and energy are limited.The density is limited to some minimum fraction and a maximum multiple of the average ghost-zone value at that level.The energy per unit mass is limited to a range of values about its boundary value.The density derivative at the boundary is further modified by replacing the vertical derivative of vertical momentum with the momentum divided by the scale height: where H P is the pressure scale height and y is the vertical direction pointed inward.The ghost-zone velocities are set to their value at the boundary and are limited to a maximum absolute value for both inflows and outflows, and which cannot be greater than their value at the boundary for inflows.
The magnetic field is made to tend toward a potential field at the boundary.This is achieved by creating a potential field extrapolation, B potential , of the magnetic field into the ghost zones and using this to set the vertical derivative of the horizontal components of the electric field proportional to ( ) B B y t potential -D D .

Input/Output
A snapshot consists of the nine variables: density, three components of momentum, internal energy per unit volume, temperature, and three magnetic field components.If a from="from.dat"file is specified in the &io namelist, it is used to read an initial snapshot specifying all nine variables that will be calculated evolving in time, otherwise the same file="run.dat"is used for both input and output.In addition, a mesh file must be selected that gives the coordinate, coordinate spacing, and Jacobian between physical and index space at each grid location.The first line gives the three dimensions.Then, for each dimension there is a line giving the following: the cell-centered grid spacing, the face-centered grid spacing, the cell-centered location, the face-centered location,  6, 630.15, 630.25, 1564.8, and 1565.2 nm lines and the Si I 1082.7 nm line with the German Vacuum Tower Telescope and Hinode (Beck et al. 2013).The continuum intensities were matched except for the 557 nm line, which is in absolute flux units.The SiI line is formed in the upper photosphere where both non-LTE effects neglected in the simulations as well as boundary effects reduce the agreement with the observations.the 1/derivative of the cell-centered locations with respect to the cell index, and the 1/derivative of the face-centered locations with respect to the cell index.If there is a from.datfile, these values will be taken from the file from.mshif it exists.Similarly, a from.timfile is sought with the time information.
The output consists of 11 files based on the root name run' of the run.datfile.One is run.datitself, to which snapshots are written every nstep time steps or tsnap time intervals.Another is a run.scrfile of snapshots that overwrite each other every nscr time steps or tscr time intervals.These are direct-access unformatted files.At the same interval as the saved snapshots are written, file run.tim has the current time information added to it, a file of the various energy fluxes is written to run.flux, and various statistical quantities are written to run.stat.A log file is written for each run.The script that starts the run begins it by recording the date and time of the start of the run, text describing the run, the size and date of the executable, the number of CPUs per node, the number of nodes, and a list of the nodes.As each namelist is read the value of each parameter, whether default or read in from the namelists, is printed to the log file.Next, if a snapshot of the variables is read the name of the "from" file is given, then the time, if a time file is found, and then the first several values of each variable near the top (lb) is printed as they are read in.At every time step the log file is appended with the time step number, the time, the time interval, and the location of the most restrictive of each of the Courant conditions due to viscosity, resistivity, density, and energy changes.At each time step a small text file run.chk is also written that gives similar time step information.At the start of each run a copy of the input file is saved as well as added to a file of all the input files for that run, a copy of the mesh file is saved, and a file giving the dimensions, grid spacing, number of variables, and MPI decomposition is written.

Tests and Comparisons with Observations
We have performed several basic tests of the Stagger code as well as comparing its results with observations of solar spectral lines, the solar velocity spectrum, and solar p-mode properties.

Advection
To test the advection properties of the Stagger code, we advected a density hat with a constant velocity toward the right, in Figure 5.The smooth density profile is barely altered after being advected 480 grid cells.The steep density profile is advected properly but develops oscillations behind the sharp fronts.The Stagger code is conservative (the total mass is unchanged) and nondissipative.The Stagger code is, however, dispersive; that is, shorter-wavelength components travel slower than the longer-wavelength components, which leads to short-wavelength oscillations in the density after advection.Pereira et al. (2013).Top panels: comparison with the visible/infrared observations of Neckel & Labs (1994).Bottom panels: comparison with the near-infrared observations of Pierce & Slaughter (1977) and Pierce et al. (1977) for wavelengths between 1158.35-2401.8nm.
The conservative, nondissipative behavior is not affected.In order to minimize the dispersive effects, it is necessary to use sixth-order derivatives and fifth-order interpolations of the density on the staggered grid.In the presence of steeper fronts (shocks) these dispersive effects appear, and extra diffusion in Equations (2), (3), and (5) is needed to damp out the short, unphysical wavelengths.

Brio-Wu MHD Shock Tube
A Brio-Wu shock tube (Brio & Wu 1988) was calculated.The initial state (Table 2) is a single step in density, pressure = energy (γ = 2) and transverse magnetic field, near the middle of the domain.The velocity is zero everywhere to begin with.
Results were computed on a grid of 250 cells and compared with a calculation using 8000 cells at time t = 0.15 (Figure 6).A fast rarefaction is the first wave to the right.It has a small overshoot in the velocities at low resolution.Behind is a shock that is well handled at both resolutions for all the variables and which moves slightly faster at low resolution.Density and momenta have jumps at the contact discontinuity, which is broader at low resolution.A slow compound wave and rarefaction farther to the left are slightly broader at low resolution.There is also some rounding at transitions at low resolution.
Convergence properties are shown in Figure 7. Comparing density profiles at different resolutions shows that resolution primarily affects the various fronts, while smooth portions of the profile are nearly resolution independent.The difference in the integral of the density and absolute value of the velocity at low resolutions with that at the highest resolution shows that the convergence at these fronts is approximately linear in the grid spacing.

Orszag-Tang
The Orszag-Tang 2D MHD test was also run (Orszag & Tang 1979).It tests how well a code handles the formation of MHD shocks and shock-shock interactions, and also reveals if any directionality is imparted on the dynamics by the grid.The initial state is as follows: The resulting gas and magnetic pressure are shown in Figure 8.The performance of the Stagger code is evaluated by comparing these results with similar calculations using other codes, e.g., Dai & Woodward (1998).
Shocks are visible as sharp jumps in the pressure (Figure 8).Interactions of shocks and waves are further illustrated in Figure 9, showing contours of the negative of the velocity divergence and of the current density superimposed on the gas pressure image.The shocks around the central oval and the long (red) divergence trails from it are fast shocks (Snow et al. 2021).

Spectral Lines
A sensitive test of both the dynamic and thermal structure of the simulations is given by comparing observed solar spectra with spectral synthesis performed on a realistic solar simulation.The spectral energy distribution (i.e., low-resolution spectrum) from 300 to 1000 nm is in good agreement with observations (Figure 10), showing that the mean thermal structure of the Stagger solar simulations is realistic.
Spectral line profiles depend sensitively on both the temperature structure and the convective velocities.The line profiles at any given location vary in depth, wavelength, and width.To compare with observations, the simulated profiles are convolved with a spatial Gaussian times a Lorentzian kernel (representing the telescope point-spread function) and convolved with a Gaussian in wavelength (representing the instrument degradation) to obtain the same resolution as the observations and then spatially averaged to produce the mean profile.These mean profiles are an accurate check on the detailed correlation of temperature and velocity in the solar Stagger simulations.Figure 11 shows comparisons with several iron lines in the visible and infrared formed in the low to midphotosphere (Beck et al. 2013).The simulation covered a region 6 × 6 Mm wide and extended from the temperature minimum down to a depth of 2 Mm below unit continuum optical depth.The horizontal grid spacing was 24 km and the vertical grid spacing was 15 km.The simulations used four (continuum, weak, medium, and strong lines) scaled opacity bins (see Section 7.7.1).The continuum and line opacities and the EOS for the simulation were taken from the Uppsala stellar atmospheres package and ODF tables (see Section 8.1).The spectral line synthesis was subsequently performed with the LTE LILIA code (Socas-Navarro 2001), based on a number of simulation snapshots.It is the convective velocities that provide the nonthermal Doppler broadening in the simulation with no free parameters, which in 1D calculations is represented by micro and macro turbulence parameters.
The line core velocities also match well between the degraded simulation and the observations (Figure 12).Note that the spectrum is featureless except for the peak at granulation scales, ∼1 Mm, and excess power at supergranulation scales.The amplitude at very low degree is uncertain.
Further comparisons between the simulation-generated and observed spectral lines, e.g., intensity histogram, line depth, line asymmetry, equivalent width and FWHM, can be found in Beck et al. (2013).With the current, larger New Solar Telescope and DKIST solar telescopes more detailed comparisons will be possible since the telescope resolution reaches that of the simulations (24 km).A few simulation snapshots with 12 km horizontal resolution are available, but may not be adequately relaxed to the higher resolution.

Center-to-limb Variation
The center-to-limb variations of continuum intensities provide a sensitive probe of the solar photosphere.Because the continuum intensity is proportional to the local source function of continuum-forming regions, its center-to-limb variation is a measure of the temperature variation with depth (the closer to the solar limb, the higher up in the atmosphere).Pereira et al. (2013) have made such a comparison, shown in Figure 13.For this work, Stagger simulations based on the EOS and opacities detailed in Sections 8.1 and 8.2 were compared, both adopting the solar chemical composition by Asplund et al. (2005).The emergent intensities were spatially and temporally averaged over 90 snapshots.There is a significant improvement between the older four-bin scaled opacities and the newer 12bin optimized opacities that are divided in both opacity strength and wavelength (Figure 2).

Convection Velocity Spectrum
The convective velocity spectrum is a good test of the velocity amplitudes at different spatial scales.This in turn gives information on the velocities at different depths, since the scale of the dominant convective structures increases with depth because of mass conservation in the stratified atmosphere.Figure 14  To compare with the observations, synthetic spectral lines were calculated for the vertical columns in the simulation results using the Nicole code (Socas-Navarro et al. 2015) assuming LTE.These are then modified to reproduce the instrumental degradation of the observations.The vertical velocity is determined from a Doppler shift of the spectral lines and the horizontal velocities are determined by local correlation tracking (LCT) applied to successive continuum images for both the simulation and the observations (Figures 15 and 16).
The velocities are then compared with the actual simulation velocities on surfaces of constant optical depth (τ 500 = const.)corresponding to the formation region of the spectral line.The best correlation for the vertical component (0.96) is obtained at τ 500 = 0.1.The best correlation for the horizontal component is lower (0.5) at τ 500 = 1, which provides an estimate of the accuracy of LCT.

p-modes
The p-mode spectrum samples the mean structure throughout the simulation domain.We restrict comparisons with solar observations to nonradial modes that have their lower turning points well within the computational box of 20 Mm depth.
The solar MHD convection simulations were run with the mconv experiment.We make a direct comparison of an 18 hr time sequence of a 96 Mm wide by 20.5 Mm deep simulation with a 36 hr 93.3 Mm square patch of the HMI Doppler velocity observations at disk center.For these simulations, we  used the four-bin scaled opacities and an EOS, both from the Uppsala stellar atmospheres package (Gustafsson et al. 1975).The simulation had 1 kG horizontal magnetic field advected into the computational domain by inflows at the bottom.At continuum optical depth 0.1, the average unsigned vertical field is 2 G and its maximum is 2 kG in either direction.The average horizontal field is 3.5 G with a maximum of 650 G.We take the Fourier transform in space and time of both data sets.Because of the slight difference in horizontal dimensions, we linearly interpolate the Fourier-transformed HMI data to the same wavenumber grid as the simulations.The k-Ω diagram is shown in Figure 17.
Individual spectra for ℓ = 410, 546, and 729 are shown in Figure 18.The good, but not perfect, agreement in mode frequency, width, and amplitude between the simulation and the observed modes indicates that the simulated convective structure is an accurate representation of the solar structure in the top 20 Mm of the convection zone and that the driving and damping processes are fairly accurately represented in the simulations.The small deviations provide important data for further analysis of p-mode driving and damping.Theoretically (Goldreich et al. 1994;Stein & Nordlund 2001), we expect the modes to be excited stochastically by rapid variations in the Reynolds stresses and surface entropy fluctuations.The mode widths represent the mode lifetimes, which range from 1 to 4 hr for these modes.Lifetimes should depend primarily on the mode damping.The significant deviations between the simulated and observed modes are due to the following.(i) The lowest-frequency simulated mode amplitudes are larger than those observed for all wavelengths.(ii) Simulated mode widths increase with decreasing wavelength relative to those observed, from slightly narrower than observed to wider than observed.(iii) The simulated mode amplitudes at each wavelength decrease relative to those observed with increasing frequency and become equal for the second or third harmonics.(iv) The intermode power in the simulations is larger than observed.

Scaling
We have tested the scaling of the Stagger code on the Blue Waters and Pleiades supercomputers using a 48 × 48 Mm wide by 20.5 Mm vertical computational box of dimensions 4032 2 × 500.It included radiative transfer with a four-bin scaled opacity, and one vertical and four slanted rays.The runs were for 100 time steps, using 4032 MPI tasks.Figure 19 shows the wall time per time step for 3200-64,000 processors.
Where there are multiple data points, they refer to multiple runs with the same setup.The difference in times is due to different loads on the computers.The difference between Pleiades and Blue Waters is due to the different processors on the machines.Significant departures from linearity occurs at 50 3 grid points per processor.

Conclusion
The Stagger code is a realistic 3D MHD code.It is an excellent tool for modeling many astrophysical situations and has been used to study stellar atmospheres, star formation, galaxy evolution, and supersonic turbulence.It is modular and can easily be modified to include the needed physics.Such codes can be used both to make predictions that can guide and be confirmed by observations (e.g., edge brightenings of solar granules were predicted before they were observed) and aid in the interpretation of observations by producing models with much more information than is observable.In this way, our understanding of astrophysical phenomena is increased.
In this paper, we have given a detailed description of the methodology used in the Stagger code: the sixth-order centered derivatives and fifth-order centered interpolations on a staggered grid; the third-order, low-memory-usage Runge-Kutta time advance; and the various methods of applying boundary conditions using ghost zones in the nonperiodic direction.The different modules for treating energy transfer were presented, including a detailed description of the nongray radiative transfer solution.The Feautrier equations are solved on a rotating set of long characteristics, with one vertical and several slanted rays.Opacities are binned to drastically reduce the number of frequencies solved for.LTE is assumed.We described the physics used in the EOS and opacity calculations.Parameters are input via a text file containing a number of namelists with their values.We provide several tests of the Stagger code, both for idealized problems and in comparison with solar observations.These show that the Stagger code is highly conservative, nondissipative, with low dispersion.The behavior of the dynamic surface layers of the solar atmosphere is accurately modeled.The Stagger code has excellent parallel scaling, including with a realistic nongray radiative transfer solution, when there are more than 50 3 cells per processor.
Much progress has been made in astrophysics with past and current computing capacity.Exabyte computing capability and new codes to take advantage of them will significantly increase the extent and complexity of astrophysical phenomena that can be simulated.
routines that read each namelist.If the parameters in a given namelist are not needed, the namelist can be left blank, or with only selected parameters present.In addition, there are some initialization routines that may or may not be compiled into the executable that have additional namelist inputs.The most important namelists are shown below, with typical input parameters commented upon.There are also other, more specialized versions of this namelist.
Similar equations are used in each dimension and for downward interpolations.As for the derivatives, MPI routines are used to get values at the opposing periodic boundaries, and ghost zones are used for boundary values in nonperiodic directions.

Figure 1 .
Figure 1.Computational cell showing the location of the physical variables.

(
in, for example, accretion disk simulations, to allow for a hot corona around the disk.It is also sometimes used if only a crude chromospheric model is desired.Then,Schrijver & Zwaan 2000)

Figure 2 .
Figure 2. Illustration of the opacity binning of a solar simulation, showing opacity strength, Equation (44), as function of wavelength for the 107,855 wavelengths included in this monochromatic calculation.Bin membership is shown with color, and the borders between the bins are also shown with dashed black lines.Various continuum opacity contributions are recognized in the lower envelope of these points, including the H − bump centered at log 3.9 l = , the free-free absorption redward of that, the molecular bands around log 4.2 l = , the hydrogen Lyman edge just short of log 3 l = , and the absorption edges of Al I and Ca I at log 3.3 l = , etc. Left panel: the four equidistant bins of Section 7.7.1.Right: a 12-bin implementation of the optimized binning of Section 7.7.2.

Figure 3 .
Figure3.A comparison of the 1D radiative cooling of the full monochromatic solution for a solar simulation, and the cooling for the binned solutions shown in Figure2.The differences between the monochromatic and the binned solutions are shown, exaggerated by a factor of 10, with the red and blue dashed lines on the right-hand scale.Left: the four equidistant bins of Section 7.7.1, also shown in the left-hand panel of Figure2.Scaled opacity calculation shown in blue, and individual bin-wise opacities, as in Section 7.7.2, in red.With scaled opacities the maximal difference with respect to the monochromatic calculation, relative to the magnitude of the cooling peak (shown with vertical dotted line), is 21.75%, the overall rms deviation is 6.58%, and the rms deviation just in the cooling peak is 10.32%.For individual opacities the deviations are 2.81%, 0.63%, and 1.13%.Right: a 12-bin implementation of the optimized binning of Section 7.7.2, also shown in the righthand panel of Figure2, giving deviations of just 1.22%, 0.26%, and 0.45%, respectively.

Figure 4 .
Figure 4. Extrapolation of the 1D mean stratification to the 2D table.The points are a random subset of the simulation cube&#39;s upper 17 layers, and the lines are the fits described in the text, with coefficients a 1 -a 5 .This example is for the simulation t47g22m+00 at t = 10 6 s.
)), ρκ B (Equation (41)), and the 5000 Å monochromatic opacity, ρκ 5000 .The current tables have a resolution of N T = 221 and N ρ = 157, with an extent of [ t bdry is a relaxation timescale for the boundary.Fluctuations in the horizontal momenta are similarly damped.

Figure 5 .
Figure 5. Density advection of a smooth profile (left) and a sharp profile (right).The solid line (and squares) is the initial density distribution displaced 480 grid cells; the diamonds are the density after being advected 480 grid cells to the right.

Figure 6 .
Figure 6.Brio-Wu shock tube results for density, pressure, parallel and transverse velocities, and transverse magnetic field at time t = 0.15 computed on a grid of 250 cells (diamonds) compared to a calculation with 8000 cells (solid lines).The viscosity parameters were c v = nu1 = 0.04, c s = nu2 = 2.5, and c w = nu3 = 0.04 (Equation (20) and the Appendix).

Figure 9 .
Figure 9. Image of the pressure with superimposed contours of the negative velocity divergence (left) and current density (right).

Figure 10 .
Figure 10.Comparison of the simulated solar flux (black) from 300 to 1000 nm with the observed (red; Kurucz 2005) solar surface flux (Chiavassa et al. 2018).The good agreement shows that the thermodynamic structure of the simulation near the surface is realistic.

Figure 11 .
Figure 11.Comparison of line spectra calculated from a hydrodynamic 6 × 6 Mm Stagger simulation and observations of the Fe I 557.6, 630.15, 630.25, 1564.8, and 1565.2 nm lines and the Si I 1082.7 nm line with the German Vacuum Tower Telescope and Hinode (Beck et al. 2013).The continuum intensities were matched except for the 557 nm line, which is in absolute flux units.The SiI line is formed in the upper photosphere where both non-LTE effects neglected in the simulations as well as boundary effects reduce the agreement with the observations.

Figure 12 .
Figure 12.Histogram of line core velocities from the simulation (red) and observations (black).

Figure 14 .
Figure 14.Simulated convective total velocity spectrum compared to the observed solar spectrum as calculated by Hathaway et al. (2015) from HMI data.Plotted is [kP(k)] 1/2 , where P(k) is the velocity power at wavenumber k.Note that the spectrum is featureless except for the peak at granulation scales, ∼1 Mm, and excess power at supergranulation scales.The amplitude at very low degree is uncertain.
compares the spectrum of simulated velocities at optical depth unity with those determined by Hathaway et al. (2015) from Helioseismic and Magnetic Imager (HMI) data.The simulation data are from a 96 × 96 Mm wide and 20 Mm deep run covering 25 hr of solar time.The horizontal resolution was 24 km and the vertical resolution varied from 12 km near optical depth unity to 70 km near 20 Mm depth.A horizontal magnetic field was advected into the computational domain by inflows at the bottom boundary.The field strength was slowly increased to 1 kG and thereafter held fixed.Thus, magnetic flux is constantly being advected into the simulation domain.There is no net vertical magnetic field at any horizontal plane in the simulation, but of course in magnetic concentrations there are local strong vertical fields.There is excellent agreement in the range of granulation.However, the simulation does not show the observed excess power at supergranule scales.The reason for this is still an issue for ongoing research.Yelles Chaouche et al. (2014) compared the simulation vertical and horizontal velocities with IMaX observations at disk center of the Fe I 525.02 nm line.The simulations had a box size of 6 × 6 Mm horizontally and extended from 0.43 Mm above the average τ 500 = 1 level to 2 Mm below this level.Runs were made with zero magnetic field and with 50, 100, and 200 G average signed vertical field.

Figure 16 .
Figure 16.Horizontal velocity spectrum computed by LCT of simulation compared to IMaX-observed spectrum (Yelles Chaouche et al. 2014).The left figure is for the hydrodynamic run and the right shows results for the three MHD runs as well.The color scheme is the same as for Figure 15.

Figure 17 .
Figure 17.k-Ω diagram of p-mode power.Images are from an 18 hr time sequence of a 96 Mm square by 20.5 Mm deep simulation and contours are from a 36 hr time sequence of 93.3 Mm square HMI data.The abscissa is ℓ and the ordinate is the frequency in millihertz.

Figure 18 .
Figure 18.The p-mode spectra comparing the 96 Mm wide simulation with the 93.3 Mm wide HMI observations for horizontal wavelengths ℓ = 410, 546, and 729 corresponding to horizontal wavelengths 11, 8, and 6 Mm respectively.The simulation was an 18 hr time sequence and the HMI observations were a 36 hr time sequence.The power is summed over all directions at each wavenumber.The units of power are m 2 s −2 .

Figure 19 .
Figure19.Wall time per time step as a function of the number of processors (in units of a thousand) for a 4032 2 × 500 solar convection simulation.

Table 2
Initial State