Do MURaM and STAGGER Simulations of Solar Faculae Match Observational Signatures from Magnetic Structures?

We compare results of simulations of solar facular-like conditions performed using the numerical codes MURaM and STAGGER. Both simulation sets have a similar setup, including the initial condition of ≈200 G vertical magnetic flux. After interpolating the output physical quantities to constant optical depth, we compare them and test them against inversion results from solar observations. From the snapshots, we compute the monochromatic continuum in the visible and infrared, and the full Stokes vector of the Fe i spectral line pair around 6301–6302 Å. We compare the predicted spectral lines (at the simulation resolution and after smearing to the HINODE SP/SOT resolution) in terms of their main parameters for the Stokes I line profiles, and of their area and amplitude asymmetry for the Stokes V profiles. The codes produce magnetoconvection with similar appearance and distribution in temperature and velocity. The results also closely match the values from recent relevant solar observations. Although the overall distribution of the magnetic field is similar in both radiation-magnetohydrodynamic (RMHD) simulation sets, a detailed analysis reveals substantial disagreement in the field orientation, which we attribute to the differing boundary conditions. The resulting differences in the synthetic spectra disappear after spatial smearing to the resolution of the observations. We conclude that the two sets of simulations provide robust models of solar faculae. Nevertheless, we also find differences that call for caution when using results from RMHD simulations to interpret solar observational data.


Introduction
The strong (kilo-Gauss) magnetic fields in the solar photosphere are organized in structures that vary significantly in size, brightness, and total magnetic flux, forming a hierarchy of solar magnetic patterns (Zwaan 1978). The small magnetic structures are typically bright in continuum images, while the larger ones appear dark. The small structures are always present in the Sun (irrespectively of its magnetic cycle phase) and typically situated between convective granules (see de Wijn et al. 2009). In the quiet Sun, they appear as individual magnetic elements while in active regions they are arranged into chains, crinkles, or small clusters named filigree by Dunn & Zirker (1973). In active regions, their collective appearance is known as solar faculae when seen in the photospheric diagnostics and as solar plage when seen in the cores of the strong chromospheric lines.
The knowledge of the solar photosphere and of its magnetic structure and evolution has been dramatically advanced by three-dimensional (3D) numerical radiation-magnetohydrodynamic (RMHD) simulations (see extensive reviews by Nordlund et al. (2009) and Stein (2012)). After the pioneering work of Nordlund (1982), numerical simulations of nearsurface convection rapidly developed in the last decades in terms of the complexity of the included physics (see, e.g., Stein & Nordlund 1998;Vögler 2003;Vögler et al. 2005). In these and similar (magneto-)convection simulations, a complex set of non-local physical equations (Navier-Stokes equations of fluid dynamics, Maxwell's equations, the equation of state, and the radiative transfer equation) is numerically solved on a geometrical grid by using physically meaningful boundary conditions and very few free parameters (e.g., effective temperature, gravity, metallicity).
The results of simulations of various solar RMHD regimes are confronted with observations in various diagnostics. Typically, very good agreement is found. Some impressive examples of the match between quiet-Sun simulations and observations are the predicted shapes of photospheric spectral lines (e.g., Asplund et al. 2000a), quiet-Sun granular root-mean-square (rms) contrast (e.g., Danilovic et al. 2008), and center-to-limb variation of the continuum intensity (e.g., Pereira et al. 2009).
The 3D RMHD numerical simulations designed to study faculae/plage are commonly initiated with a vertical unipolar hecto-Gauss field providing a magnetic flux similar to that in the active regions outside sunspots (Bercik 2002;Bercik et al. 2003;Vögler 2003;Vögler et al. 2005). Keller et al. (2004) and Carlsson et al. (2004) used snapshots from this type of simulations to explain and reproduce the appearance of bright granulation walls and faculae in solar observations toward the limb. Beeck et al. (2012) performed an important cross-validation of the results of hydrodynamical simulations of solar convection obtained with three different R(M)HD codes: CO 5 BOLD (Freytag et al. 2002;Wedemeyer et al. 2004), MURaM (Vögler 2003;Vögler et al. 2005) and STAGGER (Nordlund 1982;Stein & Nordlund 1998). They compared the vertical stratification of several variables computed by the codes as well as the center-to-limb variation of the continuum intensity at various wavelengths synthesized on the basis of the snapshots. Despite the different simulation setups used to obtain each of their data sets, they found that the results for the mentioned codes agree quite well in the variables they studied. Their results are encouraging in terms of applications of the simulation snapshots, giving additional confidence in the reliability of the numerical simulations themselves. However, Beeck et al. (2012) focused only on hydrodynamical simulations, and the same applies for the detailed comparison by Kupka (2009) of results from different stellar convection codes used to produce 3D atmospheres for stellar structure models and asteroseismology.
Including magnetic fields adds complexity to the problem. Observational evidence shows that magnetism is an essential part of solar phenomena (Trujillo Bueno et al. 2004), and recent theoretical work confirms that the solar surface and upper convection zone are inherently magnetized (Khomenko et al. 2017). When comparing non-magnetic and magnetic simulations, considerable differences have been found in the flow structure in layers up to the temperature minimum (Moll et al. 2012). Thus, especially for studies requiring highly realistic model atmospheres (e.g., for accurate determination of the solar chemical composition, or for studies of solar irradiance and its variability), accounting for the effects of solar magnetic fields is a necessary step. For example, Fabbian et al. (2010) for the first time showed that direct (magnetic) and indirect (thermal) effects in 3D RMHD simulations significantly impact the derived solar composition. In the context of solar irradiance studies, faculae are believed to give the major contribution to solar variations on timescales of solar rotation. The result of 3D RMHD solar simulations suggest that photospheric magnetic fields, depending on their spatial distribution and on the wavelength considered, can affect the solar continuum intensity and flux (Vögler 2005;Fabbian et al. 2012;Danilovic et al. 2013;Criscuoli & Uitenbroek 2014a;Fabbian & Moreno-Insertis 2015;Criscuoli et al. 2020). Along with the inclusion of all important atomic and molecular spectral lines, a description based on 3D RMHD photospheric simulations (rather than semi-empirical models), should thus further improve the understanding of solar irradiance variations on different timescales (and, thus, of solar irradiance reconstruction models), and of the variable magnetic activity of solar-type stars (see, e.g., the review by Fabbian et al. 2017). Buehler et al. (2015) obtained properties of the magnetic fields in an extended solar facular/plage region from a spatially coupled inversion of observational spectropolarimetric data from HINODE. They used spatially coupled inversions, finding that plage magnetic flux concentrations (MFCs) are surrounded by a ring of strong downflows (which shows strongly broadened Stokes profiles and often harbors magnetic patches of opposite polarity compared to its corresponding main magnetic flux concentration). They also found that the MFCs expand with height, forming magnetic canopies having weaker and more inclined magnetic fields (even more inclined in the presence of large structures such as sunspots). In their inversion results, they differentiated facular pixels into those with magnetic fields respectively belonging to the (warmer) strong almost vertical magnetic core and to the (colder) weak inclined magnetic canopy.
In the work presented here, we aim to extend the study by Beeck et al. (2012) by comparing snapshots from MURaM and STAGGER 3D RMHD simulation with a mean magnetic field of the order of 200 G, i.e., representing solar facular-like conditions. The setups of the two simulation runs selected for comparison are more similar to each other than those selected by Beeck et al. (2012), though not identical.
Although it would be interesting to study setups as close as possible, such comparison would address more the numerical side of the two codes and their particular algorithms and implementations. As far as we are aware, there is no comprehensive study on the effect of different numerical algorithms, mathematical implementations, equation of state, opacities, and boundary and initial conditions in 3D solar/ stellar atmospheric magnetoconvection simulations. While this would certainly be interesting, it remains not only outside the scope of this work, but is in fact prohibitive (even if investigating this issue for a single code) due to the extremely time-consuming task of repeating such computationally intensive 3D RMHD calculations for a systematic test of each parameter. A few major efforts have, however, been accomplished at least for the non-magnetic case, providing important indications and hinting the way forward. For an extensive discussion of different boundary conditions, model relaxation times, grid resolution and initial conditions for the hydrodynamic flow, and of their consequences, we refer the reader to, e.g., Grimm-Strele et al. (2015). Our aim here is rather to study how consistent are typical 3D MHD simulations (and results derived from them) obtained by different groups. Therefore, we opted to analyze snapshots produced by MURaM and STAGGER when each of them is used with its standard choice of free parameters and supplemental data.
In Section 2, we briefly introduce the two codes and the code used for line synthesis. The details of the simulation setup and the line synthesis are specified in Section 3. In Section 4, we compare the quantities computed by the two codes to each other and against the inversion results of Buehler et al. (2015). The synthetic spectra computed from the two codes are compared and discussed in Section 5. We compare the radiative diagnostics averaged over the horizontal extent of the snapshots, then the spatially resolved diagnostics at the simulation resolution, and finally the diagnostics after they have been smeared to mimic observations at a real telescope. Finally, in Section 6, we give the conclusions of our work.

Tools
In the following subsections, we briefly introduce the tools used in this work: namely, the numerical 3D RMHD simulation codes MURaM (Vögler 2003;Vögler et al. 2005) and STAGGER (Galsgaard & Nordlund 1996), the spectral synthesis code NICOLE (Socas-Navarro et al. 2015), and the procedure for the spectral and spatial degradation for the line profiles computed from the snapshots produced by the two codes, to the resolution of the SP/SOT on board the Hinode satellite (Lites et al. 2001;Kosugi et al. 2007;Tsuneta et al. 2008).

MURaM
The MPS/University of Chicago Radiation Magnetohydrodynamics (MURaM) code was developed to study the interaction between convective fluxes and magnetic field in the solar photosphere. The version of the code used to produce the snapshots we analyze here is the one described in Vögler (2003) and Vögler et al. (2005), not to be confused with a more modern version with numerous upgrades and additional features necessary for the large-domain simulations (see Rempel et al. 2009;Rempel 2014Rempel , 2017. The MURaM simulation run is described in Vögler et al. (2005) and it has been extensively used and extended by Keller et al. (2004), Shelyag et al. (2007), Vitas et al. (2009), Afram et al. (2011, Danilovic et al. (2013), Riethmüller & Solanki (2017), and others.
MURaM uses a fourth-order central difference spatial scheme and a fourth-order Runge-Kutta temporal scheme. As described in Vögler (2003) and Shelyag (2004), the code employs a realistic equation of state (EoS) to compute temperature and pressure (needed for computing of the radiative energy exchange) in terms of density and internal energy density that are native thermodynamical quantities in the code. The EoS takes into account the effects of partial ionization. The Anders & Grevesse (1989) solar chemical mixture was adopted (both in the EoS and for the subsequent derivation of opacities). An open boundary condition was applied at the bottom of the simulation box with an algorithm implemented to control fluctuations of the total mass in the box and of the mean radiative output (see detailed description in Vögler 2003). The simulation used in this study was carried with the top boundary closed for the mass flow and with a condition of vertical magnetic field. The radiative transfer equation in MURaM is solved assuming local thermodynamical equilibrium (LTE) and using the short-characteristics formal solver (Olson & Kunasz 1987;Kunasz & Olson 1988) with 24 rays (3 per octant). The horizontal boundary condition is periodic for all the quantities. Local interpolation of the opacity is done by bilinear formula. The approximation of the quantities along the line of sight is linear. Opacity distribution functions (ODFs) used for the opacity binning (Nordlund 1982) are computed using the ATLAS 9 code (Kurucz 1993). In Table 1, we summarized the principal features of the MURaM simulation run used.

STAGGER
The 3D radiative-magnetohydrodynamics Copenhagen STAGGER code was originally developed by Nordlund (1982). For further details, see, e.g., Galsgaard & Nordlund (1996) and Beeck et al. (2012). The code is MPI-parallelized. To solve the equations for the conservation of mass, momentum, and energy, it employs a sixth-order finite-difference Cartesian scheme in space with fifth-order interpolations and a third-order Runge-Kutta scheme for time-stepping of the magnetohydrodynamic variables. The main features of the STAGGER version employed to obtain the snapshots used here are shown in Table 1.
The EoS follows Mihalas et al. (1988) and accounts for the effects of excitation, ionization, and dissociation of the most abundant chemical elements in the solar photosphere. Their solar abundance, needed in the EoS and also for the computation of opacities, is the same as used in Table 1 of Gustafsson et al. (1975).
Opacities due to continuum (Gustafsson et al. 1975, and subsequent updates) and to spectral lines (Gustafsson et al. 2008, and subsequent updates) come from the Uppsala/ MARCS opacity package. The code uses the opacity binning method of Nordlund (1982) to group the opacity sampling (OS) data into only a few opacity bins, with the assumption of LTE adopted throughout these calculations. The version of the code used to produce the snapshots we here analyze employs four opacity bins, as done in the MURaM run. The STAGGER run is described and used in Fabbian et al. (2010Fabbian et al. ( , 2012 and Fabbian & Moreno-Insertis (2015) (see also Carlsson et al. 2004 for a description of the setup) and has also been employed in Beck et al. (2013) Criscuoli & Uitenbroek (2014b), and Beck et al. (2017).
The code uses open boundary conditions vertically, both at the bottom and at the top of the simulation box. The inflow of entropy at the bottom boundary controls the radiative output in the simulation and is set as initial input by the user. At the top boundary, extrapolation of the gradient of the potential field (to ensure charge conservation) guides the magnetic field configuration. Horizontally periodic boundary conditions are imposed for all variables on the sides of the simulation box.

NICOLE
NICOLE (Non-LTE Inversion COde using the Lorien Engine, Socas-Navarro et al. 2015) is a code to perform computations for solar and stellar atmospheres, either via forward synthesis of the emerging Stokes profiles based on available physical atmospheric stratification data or derivation of atmospheric parameters via inversion based on available spectral lines. In this work, we use only the forward module of the code to compute the two Fe I spectral lines at 6301 and 6302 Å at disk center assuming LTE. The code can synthesize the spectra from 2D and 3D atmospheres, but the radiative transfer equation is solved only in the vertical direction (the 1.5D approach). The input atmosphere has to contain at least the temperature and another thermodynamical quantity (density, pressure, or electron pressure), the magnetic field vector and the line-of-sight bulk velocity. Several opacity packages and equations of state are implemented in the code. We used the opacity package and EoS of Wittmann (1974).

Simulation Setup
We used 25 snapshots computed by each of the two RMHD codes. The two chosen data sets are representative of typically used setups for each of the two codes, therefore making them a good choice for the comparison that we carry out here and useful to corroborate the results of previous works applying these codes and using similar setups.
The STAGGER snapshots were chosen as a subset of the simulation of Fabbian et al. (2010), so as to cover a total of 36 solar minutes at regular intervals of 1.5 solar minutes. In the case of the STAGGER simulation run, a uniform vertical magnetic field of ∼225 G was introduced into an already evolved radiation-hydrodynamic simulation. The MURaM simulation run was started from a plage snapshot (saved from a run that was originally initiated with the unipolar vertical magnetic field of 200 G, Vitas et al. 2009) and the output snapshots composing the MURaM simulation set we here analyze, were saved using the same time interval as in the STAGGER case.
Each simulation covers an area of ∼6.0 × 6.0 Mm 2 horizontally (i.e., ∼15 granules). Vertically, the STAGGER snapshots reach from ∼0.5 Mm above the height corresponding to the level of optical depth unity (τ 5000 = 1, where τ 5000 is the continuum optical depth at a wavelength of 5000 Å) down to ∼2.0 Mm below it. The MURaM snapshots have a vertical domain of ∼1.4 Mm, reaching ∼0.5 Mm above τ 5000 = 1.
The original STAGGER snapshots have 252 grid points (giving a uniform horizontal grid spacing of ∼24 km) in each of the two horizontal directions, and 126 non-uniformly distributed grid points vertically (with a maximum resolution in height of 15 km at the photosphere). The MURaM snapshots have 288 grid points uniformly distributed in each horizontal direction (20.8 km per cell) and 100 grid points uniformly distributed in the vertical direction (14 km per cell).

Spectral Synthesis
The output snapshots from the RMHD simulation runs with the two codes contain the mass density of the plasma, the linear momentum vector, the total energy density (for MURaM, while for STAGGER the internal energy density is saved), and the magnetic field vector. Each of these quantities is given for every vertical grid point along the common (i.e., being the same for every pixel) geometrical height scale. The quantities are interpolated to a common optical depth scale at 5000 Å and prepared in the file format required by NICOLE. In order to have a proper representation of the solar energy transport and thus of its surface magnetoconvection features, the simulations extend vertically to deep atmospheric layers. However, for the spectral synthesis, we only use the top part of the computed box (see, e.g., explanation in Section 2.1 of Fabbian et al. 2012 and in Section 2.1 of Fabbian & Moreno-Insertis 2015), since once the correct thermodynamical stratification is achieved, the physical quantities in deeper layers do not give any influence on the resulting computed spectrum.
NICOLE is then employed to synthesize the two Fe I spectral lines at 6301 and 6302 Å. The atomic data needed for the synthesis (the central wavelength λ 0 , the excitation potentials of the lower level χ e , the weighted logarithmic oscillator strengths log (gf ) and the effective Landé factor g eff L ) for the two lines are listed in Table 2. The value we adopted for the iron abundance is log(Fe e ) = 7.50 dex ).

Data Degradation
To be able to compare the computed line profiles to the profiles observed with the SOT Hinode  we convolved the computed intensities with a point-spread function (see, e.g., Nordlund 1984). We applied the ideal part of the PSF of SOT HINODE provided by Wedemeyer-Böhm (2008) and scaled to 6300 Å. We then re-gridded the smeared profiles to mimic the pixel size of the SP/SOT detector (Lites et al. 2001). Finally, we convolved each of the profiles with the instrumental profile of the spectrograph and interpolated the convolved profiles to the wavelength sampling of the instrument. The entire procedure is described in detail in Vitas (2011).

Mean Vertical Stratifications
As the first measure of similarity between the two simulation runs, we checked the average stratification of the quantities computed in the two codes. We performed the averaging over the constant τ 5000 surfaces. In Figure 1, we show the mean (isoτ 5000 ) stratification of the temperature, total magnetic field strength, inclination angle of the magnetic field, and vertical velocity. The red curve corresponds to the MURaM results, and the blue line to the STAGGER results. Throughout this work, we use the same choice of colors to show the various results of the two codes.
The mean temperature for MURaM and STAGGER is similar over the entire optical depth range here of interest for the resulting synthetic spectrum, as seen in the upper left panel of Figure 1. The MURaM mean model is cooler (by ≈100 K) at the very top layers ( ( ) t~log 4 5000 ). The two models are virtually identical for (in particular, at continuum-forming layers). The difference with the STAGGER mean model reaches up to ≈ 250 K for the deepest layers shown in the figure (i.e., for optical depth values corresponding to ( ) tlog 1 5000 ). Therefore, the MURaM mean model has slightly larger mean temperature gradient throughout the photosphere. The colored areas in this figure panel show the range of the mean temperature values in the two simulations at each optical depth. The range is similar at all optical depths. At top photospheric layers, the MURaM run reaches significantly hotter temperatures (the highest MURaM temperature being up to ≈400 K hotter than the highest STAGGER value), as well as generally slightly cooler temperatures. However, as already seen, at such high layers the MURaM average temperature tends to be slightly below that of the STAGGER run. In the midphotosphere, the STAGGER simulation reaches slightly cooler and slightly hotter values. In the deeper part of the simulation boxes, i.e., below continuum-forming layers (i.e., in the optical depth range ( ) t > log 0 5000 ), the hottest pixels of the MURaM run are hotter than those of STAGGER (while the latter still shows the coolest pixels of the two runs).
The upper right panel in Figure 1 highlights the different behavior of 〈B〉 for optical depths above the surface, where STAGGER has larger mean magnetic field than MURaM. In the increases more rapidly, and very markedly starting for layers above , with the field finally becoming nearly vertical at the top of the optical depth range. This is due to the boundary condition in the MURaM run, which forces the magnetic field at the boundaries of its simulation box to be vertical (Vögler 2003, Chapter 3.4.1), while in the case of STAGGER no condition is imposed on the field (instead, its extrapolation is naturally guided by charge conservation, see Section 2.2) and therefore 〈γ〉 tends to a value of about 45°at the top of the simulation domain. The vertical velocity (lower right panel) is averaged separately for pixels showing upflow (dotted) and downflow (dashed) motions at τ 5000 = 1 (unit optical depth). The two simulation runs have a remarkably similar mean velocity stratification in each of the two components. The mean 1D stratifications give initial confidence that the two codes provide overall similar results in the hydrodynamical part (temperature and velocity), but they also indicate rather different magnetic field topology, likely due to the different top boundary condition for the magnetic field. In the next subsection, we show and analyze the projection of the physical quantities on the isoτ surfaces to compare the two snapshots in more detail.

Maps of the Physical Variables
Figure 2 shows maps and probability distribution histograms at the τ 5000 = 1 optical depth level for the temperature, vertical component of the magnetic field, inclination angle of the magnetic field and vertical velocity for the first saved snapshots of the MURaM and STAGGER runs. The corresponding maps for the τ 5000 = 0.1 and for the τ 5000 = 0.01 optical depth levels are available in Appendix A (Figures 11 and 12, respectively). Next to each pair of maps, we show the corresponding normalized histograms to complete the picture.
We compare these maps with the findings of Buehler et al. (2015), who performed spatially deconvolved inversions of the Hinode SP/SOT observations of faculae/plage. In their data, they identified the cores of strong magnetic concentrations and the surrounding canopy. We used their criterion for pixels in the snapshots we analyzed here, namely identifying as "core" pixels those with B > 1000 G at ( ) t = log 0 5000 and magnetic field strength decreasing with height and as "canopy" pixels those having B > 300 G at and magnetic field strength increasing with height. The contours of the magnetic features plotted over the maps in Figure 2, 11, and 12 show the spatial distribution of the cores and the canopies as just defined.
In all maps, the MURaM and STAGGER snapshots show familiar features and properties of the solar convection (Nordlund et al. 2009;Stein 2012), such as hot bright expanding upflowing plasma in granules and cool dark downflowing drafts of dense (compressed) plasma in intergranular lanes. The magnetic field is organized at meso-granular scales with strong concentrations of predominantly vertical field between the granules (Cattaneo et al. 2001;Vögler 2003). Locally, this field is organized as points and ribbons having temperature at τ 5000 = 1 higher than the temperature of the surrounding downflows, but lower than the hottest granular tops. Higher up in the atmosphere, the granulation pattern reverses (Leenaarts & Wedemeyer-Böhm 2005;Cheung et al. 2007) and, at τ 5000 = 0.1, the magnetic elements are the hottest structures in the map (while at τ 5000 = 0.01 their temperature is comparable to that of the hot component of the reverse granulation). Focusing now on the temperature at τ 5000 = 1, we see ( Figure 2, first row) that the MURaM map tends to show higher temperatures. In particular, the temperature distribution histogram for all its data points (top right panel, thick solid profile) is shifted toward higher values, and it has a flat top with a hint of two populations, as opposed to the STAGGER data. The biggest difference between the two histograms is around T = 6700 K, corresponding to the granular tops in MURaM being hotter than in STAGGER, causing the excess in the mean temperature stratification (Figure 1) and the brighter continuum (see Section 4.1). The MURaM snapshot is still significantly hotter than STAGGER at τ 5000 = 0.1, but they have nearly identical distribution at τ 5000 = 0.01, in agreement with Figure 1. The change from granulation to reverse granulation qualitatively follows the change seen in Figure 13 of Buehler et al. (2015). Quantitatively, the simulations show a somewhat wider temperature range. Although the PSF had been deconvolved in the inversion process of Buehler et al. (2015), their pixel size is still considerably larger than the grid cell of the simulations (approximately factor of 15 in area) Figure 2. Maps at τ 5000 = 1 for MURaM (left column) and STAGGER (middle column), and the corresponding probability distribution histograms (right column). From the top to the bottom, temperature, vertical magnetic field strength, magnetic field inclination angle, and vertical velocity are shown. The contours show the spatial distribution of the cores (orange in the first and the last row, black in the second and third rows) and the canopies (green in the first and the last row, white in the two middle rows), as defined in our Section 4.2 and reflecting the criterion used by Buehler et al. (2015).
leading to intrinsic smearing even with δ-function PSF. The partial histograms (normalized to the corresponding total number of pixels), i.e., those obtained respectively only for the cores and only for the canopy, as identified in the simulations, are shown in the figure as thin profiles (respectively solid and dotted) for τ 5000 = 1 in the rightmost column panels of Figure 2. At the three optical heights (for the histograms corresponding to τ 5000 = 0.1 and τ 5000 = 0.01, see Figures 11 and 12 in Appendix A, respectively), the temperature distributions in the cores show excellent match between the two codes. In the canopy, there is a match in the location of the temperature distribution peak (i.e., falling at roughly 6000 K for both MURaM and STAGGER, see top right panel of Figure 2) of the histograms at τ 5000 = 1, while the STAGGER temperature distribution histogram peak is systematically cooler than the MURaM peak higher up in the model (see top right panel in Figure 11 and in Figure 12). However, remembering that we used the same criterion as Buehler et al. (2015) and thus we found the canopy contours at ( ) t =log 2.3 5000 (i. e., τ 5000 ∼ 0.005), the parameter values shown in Figure 2 correspond to material deep underneath the canopy, not to the canopy itself, while only the parameter values shown in Figure 12 (which refers to τ 5000 = 0.01) are closer to those really existing in the canopy as defined by that same criterion.
This can be directly compared to the histograms in Figure 14 of Buehler et al. (2015): for the magnetic cores, at τ 5000 = 1, the match between their results (based on inversion of solar observations) and our results (based on RMHD simulations) is excellent, with all distributions peaking at about the same value (6300 K). With increasing height, the observations suggest slightly higher temperature in the cores than the one found in the simulations. This difference may be due to the limited sensitivity of the considered Fe lines to the physical conditions in high atmospheric layers, which, in the case of the spectral inversion technique used in Buehler et al. (2015), might have resulted in reduced accuracy of their derived atmospheric parameters.
The vertical component of the magnetic field ( Figure 2, second row) has similar distribution at τ 5000 = 1 in both snapshots. The magnetic field concentrations are dominated by the field strength of 1.5 to 2 kG. The mechanism of the field intensification had been explained by Nagata et al. (2008) as convective instability. The strongest vertical field in the STAGGER snapshot is stronger than the strongest one in MURaM, as may be expected since the field used to initiate the STAGGER run was ≈25 G stronger (see Section 3.1). From the histograms of the two snapshots, it is clear that difference in the strong-B tail is due to a small fraction of all pixels. The difference in the vertical magnetic field strength remains similar at τ 5000 = 0.1 and τ 5000 = 0.01 in accordance with Figure 1. The distribution of B VERT in the cores shows a nearly perfect match between the two codes at all three optical heights. At τ 5000 = 0.01 (close to , the actual value where the canopy was defined in Buehler et al. 2015) the two histograms show similar shape, with the STAGGER canopy having a larger relative area and slightly weaker field. Our histograms of the magnetic field in the cores can be compared to Figure 5 of Buehler et al. (2015): the range of values (between 800 G and 2.3 kG) and the median (1.5 kG) agree remarkably between the observations and the simulations.
We now discuss the histograms of magnetic field inclination, which we show in Figure 2 (rightmost panel, third row). In both snapshots, the strongest vertical field (cores, according to our definition) at τ 5000 = 1 is oriented within 25°around the vertical direction (solid thin profiles in the figure). Outside strong magnetic concentrations (canopies, according to our definition) both histograms (dotted thin profiles) reveal a varied population of inclined field. However, when looking at the total histograms (solid thick profiles), apart from their main distribution peak (similar between the two), there is a difference between the two results: for MURaM, the magnetic field inclination histogram shows a secondary distribution peak at 30°, while in STAGGER, this secondary peak is around 0°, i.e., the magnetic field is nearly horizontal for those areas. As one gradually moves up to higher layers, there is an evident difference between the two simulation snapshots: in the STAGGER snapshot, the frequency of the nearly horizontal field increases with height and clearly dominates already at τ 5000 = 0.01 with the fully developed canopy; opposite to that, in the MURaM snapshot, the inclined field becomes more inclined, but there is less of it, while the population of the field with the inclination between 50°and 90°increases significantly (compare the third row histograms in Figures 2, 11, and 12). The differences in the magnetic field inclination and its distribution are likely due to a combination of two factors: of slightly different 〈B〉 and, more significant, of the different top boundary condition for B. The latter directly affects the size of the canopy, making it easier to develop in the STAGGER case. On the other hand, the ratio of the areas occupied by the canopy and the cores is 2.52 and 2.82 for MURaM and STAGGER, respectively. It is also worth noting that MURaM at all three optical heights has more field with negative inclination (the opposite polarity). These pixels are located in the downdraft surrounding magnetic elements. The MURaM snapshot overall shows more of the small-scale structure especially at the granular edges and in between the granules, suggesting that this run is less diffusive. The distribution of the field inclination in the cores shows another good match between the two simulations at all heights. In the canopy, the field in the MURaM run is less inclined than in the STAGGER run in accordance to the already explained difference in the canopy appearance. Nevertheless, the observations of Buehler et al.  Figure 17) show a more inclined field in the cores (up to 60 • from the line of sight) and nearly horizontal field in the canopy. This is the only remarkable difference between their observation and the MHD snapshots we here considered. Note that, when comparing their results with ours, one should take into account that the inclination in their work is measured in the range from 0°to 180°.
The vertical velocities at τ 5000 = 1 in Figure 2 display the same pattern in the maps and in the histograms made for the two codes. With increasing height, the two distributions narrow down in a similar way, with the MURaM distribution being a more flat-topped histogram at τ 5000 = 0.1 and slightly shifted to more negative velocities at τ 5000 = 0.01. A comparison with the lower right panel of Figure 1 shows that this difference is mainly due to the mean MURaM downflows being faster than the STAGGER downflows. The histograms in the cores of MURaM and STAGGER match each other closely at all heights and show excellent match with the observations of Buehler et al. (2015): the magnetic cores in all three are nearly at rest while the surrounding material contains narrow downflows. These downflows are visible in our histograms in the form of a negative-velocity tail. Note that Buehler et al. (2015) use the spectroscopic convention for velocity, therefore in their case, it is the outward propagating material (upflows) that is characterized by negative (radial) velocities (decreasing distance between source and observer). Figure 3 shows the mean stratification within the magnetic cores and the canopy identified in Figure 2. All quantities are averaged over isoτ 5000 surfaces. The temperature stratifications of the cores and the canopies derived from MURaM, matches well those derived from STAGGER, with MURaM showing slightly steeper gradient in the canopy and less steep in the cores. The difference between the cores and the canopy for both codes fits well the semi-empirical models (Keller et al. 1990;Solanki & Brigljevic 1992). The strength and the inclination of the magnetic field show again significant differences between the two RMHD runs, due to different magnetic field boundary condition at the top of the respective simulation boxes. In the MURaM run, the magnetic field everywhere at the top of the domain is forced to be vertical (thus, giving small difference there between the cores and the canopy in B and γ), while in STAGGER, at the top of the domain, the magnetic field in the canopy is still significantly more inclined than the field in the cores. The vertical velocity (bottom right panel of Figure 3) for the canopy shows curves with similar shape for the two codes but with a range of values that significantly differ. The MURaM canopy has downflows that reach nearly −2 km s −1 at τ 5000 = 1, while they do not exceed −1 km s −1 in the case of the STAGGER run. For higher layers (τ 5000  0.01) in the canopy, MURaM predicts downflows and STAGGER predicts upflows, both on the order of a few hundred meters per second. The curves for the magnetic cores indicate that the velocity has smaller gradient with height in the MURaM snapshot, with a velocity difference of less than 400 m s −1 at any depth in the line formation region ( ( ) t Îlog 2, 0 5000 ).

Mean Continuum Intensity
We compute the continuum intensity in the 4000-16,000 Å wavelength range from the snapshots and, for comparison, from the Harvard Smithsonian Reference Atmosphere (HSRA, Gingerich et al. 1971). Figure 4 shows the computed continuum intensities (HSRA: black line; MURaM, red; STAGGER, blue) and (in gray color) the FTS solar atlas of Neckel & Labs (1984). Note that, as could be expected from a semi-empirical model, the HSRA intensities values are similar to those characterizing the pseudo-continuum of the FTS atlas. In the figure, MURaM closely coincides with HSRA, while the synthetic continuum computed for the STAGGER run is systematically lower than those two in the entire wavelength range, with a maximum difference of ∼10% at 4000 Å. A discrepancy of up to that order is in any case in line with the maximum deviation previously found when using a single snapshot (see Shchukina & Trujillo Bueno 2001). Moreover, the atlases themselves can be affected by uncertainties of up to a few percent (Burlov-Vasilijev et al. 1998;Fabbian & Moreno-Insertis 2015;Doerr et al. 2016).
The differences in continuum intensity may be explained by different bottom boundary conditions used in the two simulations. The mass inflow at the bottom of the MURaM simulation is regulated so that the emergent radiative flux fluctuates around a prescribed value that corresponds to the quiet Sun radiative output (6.310 10 erg s −1 cm −2 ), which explains the good match of the MURaM results with HSRA (a model representing the quiet Sun).
On the other hand, the radiative output of STAGGER is not driven by the code to the solar value, resulting instead as a consequence of the bottom entropy adopted by the user for the simulated physical processes. In this case, thus, the adoption of appropriate input bottom entropy values has to be checked a posteriori after the simulation run. In the case of the STAGGER data we employed here, the effective temperature of the resulting snapshots was within ∼±20 K of the solar value. In addition, the magnetic field adopted in these facular-like simulations is sufficiently strong to affect the granulation and partly inhibit the convective energy transfer, thus leading to a lower predicted continuum intensity than in the case of HSRA and to the output atmospheric models in the STAGGER run (having a ∼12% larger magnetic flux, see Section 3.1) generally being cooler and darker than those in the MURaM run. This explanation needs to be verified by detailed analysis of the non-local processes in both simulations, and by evaluation of the impact of differences in the opacities adopted in the a posteriori spectral synthesis as compared to those originally used by the simulation codes. This is, however, outside the scope of our present study.

Mean Stokes Line Profiles
The presence of photospheric magnetic field in the simulations affects the spectral line profiles in different ways. On one hand, the strength of spectral lines is most sensitive to temperature stratification changes (e.g., those ensuing when compared to the pure hydrodynamic case, thus with a mainly indirect effect of magnetic field on the spectral line intensity profiles). The two Fe I lines here considered have similar lower level excitation potential (χ lower ∼ 3.7 eV), several eV smaller than the ionization energy to the continuum (Fe II ground state) at ∼7.9 eV. Since iron is practically all ionized at solar photospheric temperatures, Fe I spectral absorption features arise from a minority species, so the fractional change of their strength (equivalent width) due to increasing temperature is negative (see, e.g., Gray 2008). On the other hand, some spectral lines can be significantly affected by the presence of magnetic field also in a direct manner, via the Zeeman effect. The two Fe I absorption lines here considered do have non-null Landé factors (see Table 2), so they are magnetically sensitive. However, being formed in the visible spectral range, the corresponding marginal increase of the Zeeman broadening can be neglected to first order, something that would not be possible at longer wavelengths. Regarding the absolute value of the continuum intensity around the two Fe I spectral lines we study here, magnetic field is not expected to have a significant impact as they induce negligible temperature change at visiblewavelength continuum-forming layers (Fabbian et al. 2010(Fabbian et al. , 2012Fabbian & Moreno-Insertis 2015).
The spatially averaged intensity profiles we obtained for the Fe I 6301.5 Å and 6302.5 Å lines are shown in the left panel of Figure 5, together with the observed profiles from the "Hamburg" FTS atlas of the solar disk center intensity (Neckel & Labs 1984;Doerr et al. 2016).
The wings of the STAGGER and of the MURaM Stokes I profiles of the Fe I lines show an exceptionally good match of the observed profiles. Small differences with the FTS atlas data in the line wings could be expected, since the simulations we base this study on represent a plage region with strong magnetic field and the overall kinematics somewhat differ from those of the quiet Sun, of which the FTS atlas data should be representative. The line width depends mainly on the thermal, Doppler, and Zeeman broadenings, which are related in highly nonlinear ways to the temperature structure and to the structure of the velocity and magnetic field. To disentangle the effects of these three would require extensive tests outside the scope of this investigation. The excellent agreement found seems to indicate, in any case, that the summed effects felt by the wings of the two Fe I lines due to the presence of strong magnetic field are negligible.
The line-center intensity of the synthetic Fe I profiles, instead, differs quite significantly between the STAGGER and MURaM results. The MURaM Fe I profiles show only slightly shallower absorption at line center than the FTS ones, while the STAGGER profiles are significantly weaker (almost 0.1 in normalized intensity).
The fact that the two Fe I line cores reach deeper in the MURaM case than in the STAGGER one can be understood by looking again at Figure 1. The mean temperature changes by roughly 1500 K between τ 5000 = 1 and τ 5000 = 0.01 (the approximate formation range of these Fe I lines, see Khomenko & Collados 2007). One sees, however, that the mean temperature stratification of the MURaM model has significantly higher gradient (the MURaM mean temperature reaching higher values than STAGGER in the deep layers). The cores of the two Fe I absorption features we consider sample relatively high atmospheric layers.
The slightly stronger magnetic field in the STAGGER run induces a somewhat larger split of the Zeeman components in the line profiles, making them wider and, thus, shallower. As previously discussed though, the Zeeman broadening is in itself generally small for spectral lines in the visible spectral range, and the very good match of the respective line widths in the synthetic line profiles signifies that this broadening contribution does not play a dominant role when compared to the Doppler and thermal effects. Regarding the difference between the synthetic and observed line core intensities, we first remark that we do not take into account NLTE effects (Shchukina & Trujillo Bueno 2001). Moreover, regarding the solar Fe abundance we adopted for the spectral synthesis, if we had employed the lower (by 0.05 dex) value in the latest update (Scott et al. 2015) in the series of 3D hydrodynamic solar abundance compilations initiated by Asplund et al. (2000b) and Asplund (2000), the synthetic profiles from both runs would be even weaker. Even when high-fidelity data are used i.e., with all care taken to reduce errors on the observational side, the effects of magnetic field, departures from LTE, and adopted input atomic and opacity  data make the problem of matching the solar spectrum and deriving the solar chemical composition particularly difficult.
In the right-hand panel of Figure 5, we show the Stokes V profiles of the two Fe lines computed from the two codes. The STAGGER Stokes V profiles of both spectral lines are slightly wider and stronger than the MURaM profiles due to the ∼10% stronger magnetic field in the former simulation run. We will return to this question in Section 5.5 where we discuss spatially resolved area and amplitude asymmetries of Stokes V.

Continuum Intensity Maps and Histograms
The continuum intensity maps at 6300 Å, each synthesized from one representative snapshot of the MURaM and of the STAGGER simulation run, are compared in the upper row of Figure 6. Each of the two maps is normalized to the respective mean intensity. In the lower row of Figure 6, the continuum maps are shown after smearing and re-binning to match the Hinode resolution. The color scale and the range of values are identical for all panels. The respective histograms are shown next to the maps (panels in third column of Figure 6).
At the simulation resolution, the two snapshots have a very similar range of normalized continuum intensities, with the darkest pixels at a value of 0.6 and the brightest pixels at ≈1.5 of the mean intensity. The rms contrast of the continuum at 6300 Å is 14.4% for MURaM and 14.3% for STAGGER at the simulation resolution, and 7.8% and 7.4%, respectively, after smearing. These results are close to the values found in the MURaM quiet-Sun simulations (14.4% and 7.5%, respectively) and in the Hinode observations (Danilovic et al. 2008).
The width of the histograms of the continuum intensity at 6300 Å at the full resolution is similar for the two runs although its double-peaked shape is more pronounced for STAGGER. This difference disappears at the Hinode resolution.
The brightest areas in both snapshots are small bright points located in the intergranular lanes. Comparison with the second row of Figure 2 reveals that these are the location of the strongest magnetic concentrations with B VERT up to 2.3 kG for MURaM and 2.7 kG for STAGGER at τ 5000 = 1 (2.0 and 2.2 kG at τ 5000 = 0.1, and 1.6 and 1.8 kG at τ 5000 = 0.01, respectively). In Figure 7, we compare the continuum intensity versus the total magnetic field strength at τ 5000 = 1 for the two snapshots. Both simulations show the typical hook shape (Shelyag et al. 2007;Vitas et al. 2009;Danilovic et al. 2016).
Another difference between the two snapshots at the full resolution is in the shape of the brightenings related to magnetic elements: in MURaM, these brightenings (and the magnetic elements, see Figure 2) have their peaks in the centers of the lanes, while in STAGGER, we see them often split in two ribbons with a region of reduced B in the middle. In MURaM, they also appear more wiggly, while they remain rather stretched in STAGGER. This replicates different appearance of the magnetic elements in the simulation variables and, therefore, it is likely to be due to different numerical schemes and artificial diffusivities used in the two codes. At the Hinode resolution the two snapshots look indistinguishably similar and their histograms are almost identical.

Stokes I
In Figure 8, we present maps and histograms of line-center intensity, equivalent width, and FWHM for MURaM and STAGGER simulation runs for the Fe I 6301.5 Å line. The corresponding figure for the Fe I 6302.5 Å line is shown in Appendix B (see Figure 13). To more easily compare the results, in Table 3 we list the mean value and the standard deviation of the three line parameters for each line.
The mean STAGGER line profiles (Table 3) are less deep and, therefore, the STAGGER snapshot appears brighter in the line core. The difference is especially clear in the granular tops (where MURaM is considerably darker) and in the magnetic elements (where STAGGER shows more bright points). These differences, due to the temperature gradients in the MURaM snapshot, are present in both the mean model ( Figure 1) and in the temperature maps at τ 5000 = 1 and τ 5000 = 0.1 (cf. Figures 2  and 11, note that MURaM is hotter at τ 5000 = 1 and cooler at τ 5000 = 0.1). The standard deviation of I core in both snapshots is similar.
Finally, the mean values of EW and the full-width half maximum for the two lines and the two simulation runs have same mean values and dispersion until the second decimal (see Table 3).

Stokes V
The Stokes V line profiles computed from the simulation snapshots show many deviations from the antisymmetric twolobes shape reflecting the complex variation of the magnetic field, the velocity and the temperature along the line of sight.

Number of lobes
The first indicator of the profile complexity is the number of their lobes. Khomenko et al. (2005) defined all profiles with two lobes as regular and used only those in the asymmetry analysis. The magnetic field of the individual column-atmospheres in the snapshots has many small-scale variations with height, especially at the weak-field locations, causing tiny fluctuations in the circular polarization in the essentially noisefree synthetic profiles. It is, thus, impractical to count the lobes of these profiles simply by counting the zeros of the Stokes V line profiles. Instead, prior to counting the lobes, we introduce an artificial amplitude-depending threshold to discard the small Stokes V profile variations around zero circular polarization. All Stokes V profiles are normalized to their own continuum intensity and the threshold t is set to eliminate all wavelength  0.15 ± 0.02 0.13 ± 0.03 0.15 ± 0.02 0.14 ± 0.03 points in Stokes V profiles with values t times smaller than the maximum normalized amplitude of the profile. Alternatively, we could model realistic noise with a constant level for the entire sample. However, by doing that, some of the weak Stokes V profiles and their lobes would be discarded and our aim is to compare the profiles of the two simulations also at the weak-field locations. The percentage of Stokes V profiles with a certain number of lobes at the simulation resolution is shown in Table 4. The results obtained from the two simulations are very similar. The slightly higher percentage of two-lobed, i.e., "regular" profiles in the STAGGER simulation is possibly due to slightly stronger field in that simulation. The Stokes V profiles with more than two lobes are located within the granules (the least magnetic regions) and have weak absolute amplitudes. After applying the PSF of HINODE SP/ SOT, the fraction of regular profiles in both simulations is over 99% and it remains at that level after resampling to the pixel size of SP/SOT and after applying the instrumental profile.

Stokes V asymmetries
The area (δ A ) and amplitude (δ a ) asymmetry of Stokes V profiles are important quantities indicating gradients of magnetic and velocity field along the line of sight (Illing et al. 1975;Solanki 1989). The two quantities are defined as: where a b and a r are the amplitudes of the blue and the red lobe of Stokes V profile, and V is the profile itself. Various types of Stokes V asymmetries are observed in the quiet Sun and active regions. Their complexity is interpreted in terms of one-dimensional, one-or two-component models of magnetic concentrations (Grossmann-Doerth et al. 1988, 1989Martínez Pillet et al. 1997;Bellot Rubio et al. 2000;Grossmann-Doerth et al. 2000;Steiner 2000;López Ariste 2002;Sigwarth 2001;Socas-Navarro et al. 2004) or in terms of 2D (Sigwarth et al. 1999) and 3D simulations (Khomenko et al. 2005;Shelyag et al. 2007) of the quiet Sun and solar faculae/plage. Rezaei et al. (2007), with the SP/SOT Hinode, observed the internal structure of a magnetic element and the related asymmetries, and identified the same structure in 3D numerical simulations. Using the same instrument,  observed the Stokes V asymmetries in the quiet Sun and interpreted them in terms of microstructured atmospheres. In a sequel paper,  performed the most detailed classification of the Stokes V shapes finding that 93% of all the profiles in their sample show asymmetries.
In Figure 9, we show the amplitude asymmetries (δ a ) and, in Figure 10, the area asymmetries (δ A ), as they are defined by Equations (1) and (2) for the two snapshots and the Fe I 6301.5 Å spectral line. The upper row shows the quantities evaluated at the full resolution of the simulation; the lower row the quantities evaluated after the synthetic profiles have been degraded to the resolution of SP/SOT Hinode. The corresponding figures for Fe I 6302 Å are similar (see Figures 14 and  15 in Appendix B).
Some of the profiles are so distorted from the antisymmetric shape that the above definition of the asymmetries makes little sense. A clear example of such profiles are those with two negative lobes or "two-humped profiles" (Steiner 2000); in the MURaM and STAGGER snapshots we based our study on, such heavily distorted profiles always correspond to the low magnetic field. If we nevertheless want to compare our results with those of Shelyag et al. (2007), we must directly apply Equation (1) to all synthetic profiles including the extremely distorted ones. The distribution of our MURaM results (upper left panels of Figures 10 and 15) show the same pattern as in their Figure 3 (the middle row, second and third column).
The area asymmetry is dominantly positive with small negative patches located mostly at the granular edges. On the other hand, at the simulation resolution, there is a striking difference between MURaM and STAGGER in both asymmetries and in both lines: the STAGGER profiles (the top right panels of Figures 9, 10, 14, and 15) show large patches of negative asymmetries that are entirely missing in the MURaM data. This difference is clearly visible in the corresponding histograms as well as in the mean values listed in Table 5. This table collects the mean values of the asymmetries for the cores of the magnetic field concentrations and the surrounding canopy sampled as in Buehler et al. (2015) and plotted over the maps in our Figure 2.
In the cores, the two snapshots show good match in both the area and the amplitude asymmetry, and at both resolutions. The 6301 line has larger δ a , but smaller δ A than the 6302 line at full resolution, while both lines give the same result at the resolution of Hinode (δ A ≈ − 0.01, δ a ≈ 0.1). The amplitude asymmetry in both lines is always positive. The mean area asymmetry can be both positive and negative. It has negative contribution from the magnetic core and positive from the canopy. This applies to both snapshots and in both Fe lines and it agrees with the observations of Martínez Pillet et al. (1997).
At the Hinode resolution (the lower panels of Figures 9 and  14), however, the difference in the Stokes V asymmetries between the two codes disappear: both snapshots show only small patches of the negative asymmetry in dominantly positive distribution and the corresponding histograms match excellently.
To explain the change in the distribution of the asymmetries with smearing, we added a contour to the panels with the simulation resolution showing where the amplitude of the Stokes V profiles is larger than a threshold (8 × 10 −2 ). The contours surround all the large patches of negative asymmetries that correspond to the locations of granular tops. Due to their low V-amplitudes, these profiles contribute insignificantly to the asymmetry of the profiles at the Hinode resolution. Histograms of the asymmetries at the full resolution for the profiles with the V-amplitudes above the threshold are added to the histograms in Figures 9, 10, 14, and 15 (thin curves): they show a close match between the results of the two codes.

Conclusions
We analyzed results of 3D RMHD simulations of solar faculae performed using two workhorse astrophysical numerical codes. We aim to compare the results of this type of simulation for solar facular-like conditions for the first time. Thus, we chose snapshots from simulations initiated with similar setups. In addition, we compared radiative diagnostics synthesized from these results at the grid resolution of the simulations and after spatial and spectral degradation to the resolution of the SP/SOT Hinode.
The results of the two codes show an overall agreement, confirming the conclusions of the cross-validation study that Beeck et al. (2012) performed in the case of hydrodynamical quiet-Sun simulations. Both codes show similar convective pattern with kilo-Gauss magnetic field concentrations in the intergranular space. The basic quantities (temperature, magnetic field, and velocity field) averaged over the isoτ surfaces in both codes are similar as well. The distribution of physical quantities for the samples of 1D columns (extracted from the 3D snapshots) that are representative of the cores of the magnetic field concentrations and of the surrounding canopy, is nearly identical in the two codes, as well as exceptionally similar to that derived from spatially coupled inversions of a plage region recently observed with SP/SOT by Buehler et al. (2015). In addition to the comparison that we have done here, it would also be interesting in the future to compare in detail our synthetic line profiles with those in the raw data acquired by Buehler et al. (2015) as well as those they obtained after inversion/deconvolution.
Beside the overall similarity, there are differences in the details of the results of the two simulation runs that call for caution. Although the simulations differ only slightly in numerical diffusivities, equation of state, opacities, radiative transfer solution and geometrical grid, there are significant differences in some of the boundary conditions (namely, in the way of adjusting the entropy at the bottom boundary and in the condition for the magnetic field at the top boundary), which are likely to cause differences in the results: the MURaM run is somewhat hotter across almost the whole range of optical depths considered, while the STAGGER run has a more inclined magnetic field in the mid-photosphere. These differences leave Figure 9. Stokes V amplitude asymmetry of Fe I 6301.5 Å from MURaM (left panels) and STAGGER (middle panels) snapshots at the simulation grid (upper row) and after degradation to the SP/SOT Hinode spatial and spectral resolution (lower row). The black contour separates locations with the absolute Stokes V above and below a threshold (see the text; the ticks at the contour indicate the low-amplitude side of it). Histograms in the right column compare the distribution of the amplitude asymmetry between the two codes. The thick curves count all profiles, the thin curves only the high-amplitude profiles. All histograms are normalized to the total number of pixels in each of the snapshots.
imprints on the continuum and Stokes spectral line profiles computed from the snapshots. However, after smearing the radiative diagnostics to replicate the SP/SOT Hinode observations, the differences between the two simulation runs entirely disappear. Comparison with observations obtained with the new 4 m class of solar telescopes (DKIST (Tritschler et al. 2016), EST (Collados et al. 2013)) and eventually with a nextgeneration space mission for solar data at ultrahigh spatial Figure 10. Stokes V area asymmetry of Fe I 6301.5 Å from MURaM (left panels) and STAGGER (middle panels) snapshots at the simulation grid (upper row) and after degradation to the SP/SOT Hinode spatial and spectral resolution (lower row). Histograms in the right column compare the distribution of the area asymmetry between the two codes. All histograms are normalized to the total number of pixels in each of the snapshots. resolution (Collet et al. 2016) should clarify which code approach (if any of the two) provides the most realistic simulations of solar faculae. Finally, key diagnostics that could be analyzed in future research efforts comparing the results of different 3D RMHD codes between each other and with observations include the center-to-limb variation of the emergent radiation in the continuum and in spectral lines, the linear polarization signal, and the line profiles of spectral lines that have an especially high sensitivity to certain physical parameters of the atmosphere. Figure 12. Maps at τ 5000 = 0.01 for MURaM (left column) and STAGGER (middle column), and the corresponding histograms (right column). From the top to the bottom, temperature, vertical magnetic field strength, magnetic field inclination angle, and vertical velocity are shown. The contours show the spatial distribution of the cores (orange in the first and the last row, black in the second and the third) and the canopies (green in the first and the last row, white in the two middle rows), as defined in our Section 4.2 and reflecting the criterion used by Buehler et al. (2015). Figure 14. Stokes V amplitude asymmetry of Fe I 6302.5 Å from MURaM (left panels) and STAGGER (middle panels) snapshots at the simulation grid (upper row) and after degradation to the SP/SOT Hinode spatial and spectral resolution (lower row). The black contour separates locations with the absolute Stokes V above and below a threshold (see the text, the ticks at the contour indicate the low-amplitude side of it). Histograms in the right column compare the distribution of the amplitude asymmetry between the two codes. The thick curves count all profiles, the thin curves only the high-amplitude profiles. All histograms are normalized to the total number of pixels in each of the snapshots.