Accurate Short-Characteristics Radiative Transfer in A Numerical Tool for Astrophysical RESearch (ANTARES)

We aim to improve the accuracy of radiative energy transport in three-dimensional radiation hydrodynamical simulations in ANTARES (A Numerical Tool for Astrophysical RESearch). We implement in the ANTARES short-characteristics numerical schemes a modification of the Bézier interpolant solver. This method yields a smoother surface structure in simulations of solar convection and reduces the artifacts appearing due to the limited number of rays along which the integration is done. Reducing such artifacts leads to increased stability of the code. We show that our new implementation achieves a better agreement of the temperature structure and its gradient with a semi-empirical model derived from observations, as well as of synthetic spectral-line profiles with the observed solar spectrum.


Introduction
Radiation hydrodynamical (HD) simulations have found widespread use in stellar physics in general and in solar physics in particular. Examples of their applications include the interpretation of detailed profiles of spectral absorption lines that reveal information on stellar abundances, photospheric pressure and temperature stratification, as well as convective velocity fields (see, e.g., Dravins, Lindegren, and Nordlund, 1981;Asplund et al., 2000;Nordlund, Stein, and Asplund, 2009) the characterization of the entropy jump that determines the depth of stellar envelopes and thus stellar radii (see, e.g., Ludwig, Freytag, and Steffen, 1999;Trampedach et al., 2014;Magic, Weiss, and Asplund, 2015;or Tanner, Basu, and Demarque, 2016) and also influences p-mode frequencies (near-surface effect; see Rosenthal et al., 1999;Houdek and Dupret, 2015), or, in the case of simulations that also account for a magnetic field (MHD), the dynamics of magnetic structures (e.g. Beeck et al., 2015). It is a major achievement of numerical simulations to agree in many aspects with solar and stellar observations. Since its first detailed description by Muthsam et al. (2010), the code "A Numerical Tool for Astrophysical RESearch" (ANTARES) has also been applied to a number of problems in solar and stellar physics and, with post-processing of simulation time-series to simulate observations, has allowed for much more detailed and realistic interpretation of actual observations. ANTARES has been used in studies of solar granulation statistics (Lemmerer et al., 2017), which included comparisons with lifetime and size distributions of solar granules observed with the Sunrise Imaging Magnetograph Experiment (IMaX) experiment (Martínez Pillet et al., 2011). Another application has been the analysis and physical interpretation of the damping rates of solar p-modes (Belkacem et al., 2019). The numerical simulation presented by the latter has also been used to demonstrate that patching a horizontally averaged structure based on an ANTARES three-dimensional HD model on a solar structure model is equally successful in explaining the structural part of the near-surface effect on solar p-modes as the earlier study of Rosenthal et al. (1999). Details for this can be found on a webpage about the TOP tool, which is a Python library for asteroseismology published by D. Reese, J. Ballot, and B. Putigny at top-devel.github.io/top/index.html and where examples based on ANTARES simulations are given at top-devel.github.io/top/examples.html#surface.
In general terms, ANTARES solves different variants of the equations of hydrodynamics. Among others, the code can be used to compute HD simulations of stellar surface layers, which in many cases contain convective zones situated towards the interior of the star and adjacent, locally stably stratified, photospheric layers. To this end the fully compressible Navier-Stokes equations and the stationary limit of the equation of radiative transfer are solved simultaneously and realistic microphysics is taken into account by interpolations in tables of the equation of state and opacities. Solutions for domains in one (1D), two (2D), and three (3D) spatial dimensions can be computed. If the geometrical depth of the simulation domain (spatial coordinate aligned with the radius of the star) is small enough, the spherical curvature of the star can be neglected and a Cartesian geometry can be used to define the simulation domain, which then is called a "box-in-a-star" type configuration. When coupled to solutions of the hydrodynamical equations in one, two, or three spatial dimensions, the radiative-transfer equation is always solved having a 3D configuration in mind (thus, the extra dimensions are padded with copies of the 1D or 2D models). Numerical simulations of astrophysical objects can also be performed for different grid geometries. For a general overview of the code we refer again to Muthsam et al. (2010). To perform realistic simulations of the upper part of solar or stellar envelope convection zones and their lower and middle photosphere regions, the radiative-transfer equation has to be solved. For that, a sufficiently accurate solution is needed to perform a-posteriori high-precision spectral line synthesis.
Since, for a realistic HD/MHD simulation, the most time-consuming part is the computation of the radiative transport of energy, it is worth finding the best balance between the computational speed and accuracy. In most available 3D HD/MHD codes, such as ANTARES, MURam (the Max-Planck-Institute for Aeronomy/University of Chicago Radiation Magneto-hydrodynamics code: Vögler et al., 2005), and, optionally, in CO 5 BOLD (COnservative COde for the COmputation of COmpressible COnvection in a BOx of L Dimensions: Freytag et al., 2012), a short-characteristics scheme is implemented as developed by Kunasz and Auer (1988). To save computation time, an approach of linear interpolation of the source function and opacity are used in different codes, such as STAGGER with long characteristics scheme (Magic et al., 2013a), and MURaM and CO 5 BOLD with shortcharacteristics scheme. However, according to Kunasz and Auer (1988), the computation of emergent intensity using a short-characteristics scheme with linear interpolation leads to an artificial broadening of rays. To this end, the short-characteristics method by Kunasz and Auer (1988) has originally been implemented into ANTARES with parabolic interpolation of the source function and linear interpolation of opacity. However, this scheme introduced unphysical oscillations into the solution, and ad-hoc procedures were put in place to alleviate the problems. In the mean time, more advanced variants of this method have been proposed where oscillations due to non-monotonic interpolation were automatically avoided by using the harmonic mean of left-and right-sided linear derivatives (see Auer, 2003). A modified version of this method was implemented in BIFROST (Hayek et al., 2010), however ignoring the main advantage of the method and calculating central derivatives (see Equation 15 in Hayek et al. (2010)) instead of harmonic ones of the function with height. The latter implementation still requires ad-hoc procedures to avoid unphysical oscillations, and therefore it is slower and less accurate.
In this article, we describe our modifications to improve the accuracy of the solution of the radiative-transfer equation in ANTARES and test our implementation by computing high-resolution synthetic spectra from snapshots of numerical-simulation runs, which we compare to observed solar spectra. In the following, we first describe (Section 2) how the solution of the radiative-transfer equation has been modified by us to implement an integration of the source function and of the optical depth based on quadratic Bézier splines. We then present the test of the implemented integration scheme in Section 3 and the results of our implementation (Section 4). Then we proceed with our comparison of synthetic and observed spectral-line profiles for the Sun (Section 5). At the end of the article we provide a short discussion (Section 6) and a summary (Section 7) of our results.

Radiative Transfer in ANTARES
To compute the emergent specific intensity from a 3D simulation box in which a source function and opacity are fully known, we select a number of rays that cross each grid cell and solve the radiative-transfer equation (RTE) along these rays. This numerical scheme allows the reduction of 3D radiative transfer to a set of 1D problems along such rays. ANTARES explicitly includes both upwards and downwards rays in its integration scheme.
The 1D unpolarized radiative-transfer equation in the observer's frame is where I ν is the emergent intensity, and S ν is the source function, which in local thermodynamical equilibrium (LTE) is defined by the Planck function at a specific frequency [ν]. The optical depth [τ ν ] of a path element at each frequency along a ray [r] is defined as follows: where dr = d|r | is the path length along the ray in the direction of light propagation, η ν is the total opacity at a particular frequency. The opacity is calculated as the product of an absorption coefficient [κ ν ] at a particular frequency and a density ρ [η ν = κ ν ρ]. The formal solution of Equation 1 on a discrete grid for each ray is as follows: where the index k characterizes the grid point and increases towards the inner boundary of a simulation box. Note that to simplify the notation we omit the frequency dependence in Equation 3 and much of the following, but it is still implied. Knowing the radiation field and integrating it over all rays, the mean intensity is where dω is the solid angle around the straight ray [r] along which the integration is done. Various schemes for angular integration are implemented in ANTARES (see Muthsam et al., 2010 for more details), and here we use the Gauss-Radau integration, as mentioned in Section 4.1.
Another quantity that can be calculated from the known radiation field is the radiativeenergy flux: Finally, the radiative-heating rate [Q rad ] can be computed either as or with the equivalent expression as Equation 7 follows from Equation 6 and the first angular moment of the RTE (see Sections 76 -79 in Mihalas and Mihalas (1984) or Equation 3 in Unno and Spiegel (1966)). We prefer to use Equation 6 as it avoids evaluating additional derivatives. Q rad is included in the equation of total energy conservation, which is solved alongside the other Navier-Stokes equations in ANTARES to describe the dynamics of granulation and convection at the surface of stars.

Radiative Transfer Solvers
The formal solution of the RTE (Equation 3) requires some numerical scheme to calculate the definite integral that it contains. There are different approximations for evaluating the formal solution of the RTE that are more or less accurate depending on the density of the discrete depth grid. In this subsection, we describe the RTE solvers implemented in the ANTARES code for a 1D integration along the rays. The solvers are based on a shortcharacteristics scheme (Mihalas, Auer, and Mihalas, 1978), which is defined as the line Figure 1 Radiative-transfer scheme implemented in ANTARES. The arrow shows the ray direction along which the integration is carried out in the particular case chosen here as example.
joining the intersection of the ray with the k + 1 layer, the point U (denoting "upstream"), and the current point P described by the index k (see Figure 1). Light propagates along the ray through the current grid point P . In order to determine the emergent intensity [I ] at the point P with I k ≡ I (P ), the following quantities should be specified: the intensity [I ] at the point U with I k+1 ≡ I (U), the source function and the optical depth along the ray from U to P . The emergent intensity at the previous layer I k+1 is known only at the nodes of our grid. To calculate the intensity in the intersected point U , we interpolate between the four points U 1, U 2, U 3, and U 4 using a weighted-parabolas interpolation scheme (for details see Muthsam et al., 2010). The source function and the optical depth are known only at the grid points U , P , and D (denoting "downstream") and we are free to define the interpolating functions of these quantities between U and P . We first analyze choices for the source function, in Section 2.1.1 and 2.1.2, while the optical-depth calculation is described in Section 2.2. To estimate the source function, a parabolic interpolation between U and P was originally implemented in ANTARES. To avoid complex interpolation and storage operations within the context of a 3D simulation, we require our interpolant to only need three vertical points to define it. One such scheme is a modified quadratic Bézier spline, implemented here for the source function.

Parabola-Based Solver
Originally in ANTARES, a solver using a quadrature rule according to Kunasz and Auer (1988) was implemented to integrate the source function. In this case the emergent intensity is computed as follows: where α, β, and γ are coefficients arising from an analytical integration, from k + 1 to k, of the quadratic interpolating function, going through Se −τ at the points U , P , and D. So, to obtain the intensity at P , the source function and the optical depth at the points U , P , and D should be known. The quadratic law is already more accurate compared to a linear law. However, it becomes unstable and overshoots in the presence of nonlinear gradients. In particular, this can lead to negative values for quantities that need to be positive definite, such as the optical depth [τ ] and the source function [S]. To avoid overshooting and to make the code more stable, inspections for negative values of the quantities at all points as well as I k+1 are done. In the event that a negative value is detected, it is set to zero. This numerical provision keeps the code more stable, but it reduces the accuracy of the solution of the radiative-transfer equation. In the following subsection we present our new implementation, which does not require an additional checking for any negative values and artificially setting values to zero.

The Bézier Solver
The Bézier solver is an accurate and efficient method to solve the RTE, which was proposed by Auer (2003) and implemented for unpolarized and polarized light by de la Cruz Rodríguez and Piskunov (2013). In addition to cubic Bézier spline interpolation presented by Auer (2003), de la Cruz Rodríguez and Piskunov (2013) also implemented a quadratic version of Bézier splines but showed that the accuracy of the RTE integration achieved by the two approximations is almost identical when averaging left-and right-sided control points for the quadratic spline (see below). Moreover, Bézier methods were shown to be more accurate than other methods, in particular on a coarse grid (see the analysis of Janett et al., 2017). More details on the algorithms and implementation are given by de la Cruz Rodríguez and Piskunov (2013). Because the quadratic Bézier interpolation is slightly faster than a cubic interpolation and it provides almost the same accuracy in actual atmosphere calculations (de la Cruz Rodríguez and Piskunov, 2013;Janett et al., 2017), we implement the former in the ANTARES integration scheme. With quadratic Bézier interpolation, the emergent intensity is calculated as where α , β , and γ are computed according to de la Cruz Rodríguez and Piskunov (2013).
The term with a control point C was introduced to ensure the monotonicity of the integrated function. It is computed as are the control points associated with U and P , respectively, and are defined by the function and its derivative at the appropriate grid points For the central grid point P with index k we calculate the harmonic (Bézier) derivative as Computation of the harmonic derivative leads to suppressing overshooting beyond the function values at the nodes (Fritsch and Butland, 1984). Note that for the point U the harmonic derivative cannot be computed because the function is known only at two grid points, since we want to restrict a solution at layer k to only need layers k − 1 and k + 1 for code efficiency. Therefore, we compute it from the slope between U and P . Using the average of the two control points resulting in a more symmetric spline improves the accuracy and stability of the solution. For evaluating integrals, the average of the two control points allows reaching the same order of accuracy as if a cubic Bézier spline had been used. That would not work if only one control point were used, so where possible we use two such points. But even for plain interpolation with Bézier splines accounting for the derivative of the interpolated function on both sides allows for a smaller numerical error. The stability in turn is ensured by constructing the control points such that monotonicity is ensured. Overall, this improves the stability and accuracy over methods that are either based on linear interpolations or on parabolas with ad-hoc cut-offs to suppress spurious wiggles which may lead to strictly positive quantities becoming non-positive. This modified version of Bézier integration, using the average of the two control points, we call "BézierCP2". A similar approach was applied to synthesize 3D spectral-line profiles in the IRIS-code (Ibgui et al., 2013).

Optical Depth Calculation
The numerical integration of the optical depth (Equation 2) in ANTARES was originally implemented using the trapezoidal rule. As the main goal of this article is to improve accuracy of the radiative transfer in ANTARES, we apply the more accurate Bézier integration scheme to the optical-depth calculation, too. The optical depth is written then as follows: Calculation of the control point C here slightly differs from the control points from Equation 10 and 11. Unlike the intensity integration where the emergent intensity is calculated only for the central point P , the optical depth at three grid points (U , P , and D) has to be known to solve the radiative-transfer equation. Because only differences in optical depths are required, we can always set τ at D to be zero (τ k−1 = 0). To fully implement the Bézier scheme as described by de la Cruz Rodríguez and Piskunov (2013), the opacity should be known at five grid points. This, in turn, makes the scheme computationally much more expensive. Here we implement a simplified version of Bézier integration, which uses just one control point ("BézierCP1") and not the average of two. So, the harmonic derivative is calculated only for the central grid point (P ) and then one control point is calculated for each interval (UP and P D) as follows:

Extrapolation at Grid Cell Boundaries
The solution of the RTE in the computational domain requires specified boundary conditions. The implemented boundaries are standard for radiative-transfer problems in a stellar atmosphere, i.e. no incoming radiation at the top of the domain, and at its bottom the radiation is defined by the diffusion approximation. The horizontal boundaries are periodic for all thermodynamical quantities. For more details about the boundary conditions we refer the reader to Muthsam et al. (2010) and the updates described in .
Here, we focus more on the integration at the boundaries where we can either reduce the order of the numerical scheme or extrapolate quantities to the ghost cells to implement a higher-order scheme. The latter option is used in ANTARES.
Originally, linear extrapolation of thermodynamical quantities, i.e. density and temperature, opacity, and source function to the ghost cells at the upper and lower boundaries was performed in ANTARES. As we use the diffusion approximation at the bottom boundaries, the linear extrapolation of the quantities is fine. However, at the top boundaries this extrapolation could create some artifacts and lead to an inaccurate solution. Therefore, we change the linear extrapolation of the thermodynamic quantities at the top boundary to decay exponentially for density, and we force the temperature to be constant. In instances for which the density at the lower layer is smaller than at the upper layer (ρ k+1 < ρ k ), in order to avoid exponential increase of the density at the ghost cell, we set the density at the ghost cell to be the same as ρ k . This modification improved the stability and the accuracy of ANTARES.

Test of Integration Scheme
We now compare our adopted schemes with other techniques using the example of opticaldepth integration in a 1D model atmosphere. The model atmosphere with the highest grid resolution for this test is computed using the MPS-ATLAS code for a grid with 2000 points distributed over a Rosseland optical depth [τ Ross ] varying from 10 −6.5 to 10 2 with equidistant steps in a logarithmic scale. The τ Ross , together with the opacity defines the height scale. We have computed a reference model atmosphere, and then the opacity for this model, with 2000 points per atmosphere. Then we generate other models from the reference one by taking every 20th, 50th, and 100th point. To investigate the accuracy of the integration schemes for different grids, we compare the optical depth computed in the generated models using the trapezoidal rule (called "Linear" for this test), quadratic Bézier (called "Bézier") (de la Cruz Rodríguez and Piskunov, 2013), "BézierCP1", and "BézierCP2" implemented here and described in Section 2.2 and 2.1.2, respectively. We select the optical depth computed using the Bézier scheme on the finest grid, i.e. 2000 depth points, as a reference. To compute the residuals we subtract the optical depth calculated with different methods and different resolutions from the corresponding reference value, and we normalize the results by dividing it with the reference optical depth itself and transform it into percentage values. We show the thus-normalized residuals in Figure 2. In order to capture some of the variations in a temperature gradient, we study how the optical depth and its residuals behave in the continuum and line core, which form at different heights where different temperature gradients are present. For this test we select the spectral line Fe I 6173.34 Å, which is used for the measurements of solar oscillations and magnetic field at the near surface layers (around τ ≈ 1) with the Helioseismic and Magnetic Imager (HMI: Schou et al., 2012) onboard the Solar Dynamics Observatory (SDO) spacecraft. Figure 2 shows the computed reference optical depth in the atmosphere at the continuum and at the spectral-line core together with its residuals from other methods and resolutions. As seen in Figure 2, there are two regions where the residuals increase: one region is slightly below the surface (h ≤ 0 km) and another region is at ≥ 800 km. The layers above 800 km are optically thin with τ < 10 −6 for the continuum and τ < 10 −3 for the spectral-line core (see Figure 2a) and therefore are not so important for the intensity computation. In the nearsurface region, the residuals at continuum are larger than at the line core. The reason for this is that the continuum forms in the near-surface layers where the temperature gradient varies a lot (see Figure 10) and the line core forms higher in the atmosphere and has a small contribution from the near-surface layers. For all cases in the atmospheric layers below ≈ 800 km, the "Bézier" and "BézierCP1" schemes demonstrate superior performance at spatial resolution considered. The "BézierCP2" scheme illustrates slightly larger deviation but still smaller than the linear scheme. The linear integration scheme shows the largest deviation in the entire atmosphere at all resolutions, except for the -not that important -layers at the very top where all methods converge more slowly. For the optical-depth integration in the ANTARES simulation runs, we use the "BézierCP1" numerical scheme. For additional comparisons of "Bézier" with other numerical schemes including the "Linear" one, we refer the reader to de la Cruz Rodríguez and Piskunov (2013), where the superior performance of the Bézier spline method even at moderate-to-low spatial resolution was demonstrated.
In order to test the convergence of different numerical schemes we present the root-meansquare error for the optical depth in the continuum and the spectral-line core with respect to different depth point grids in Figure 3. The larger the number of points in an atmosphere, the smaller is the error. This trend indicates that all schemes converge when we increase the resolution.

Results from Our Implementations
In this section we present the setup that we used and the main results of our study concerning the optimal radiative-transfer solution in 3D HD simulations.

Simulation Setup
We choose a cartesian box 3.68 Mm in vertical size and 6.00 Mm in horizontal size as the simulation domain. The resolution in the 2D case is 8.25 km in the vertical direction and 17.65 km in the horizontal direction. The 3D box was extended at the top in the vertical direction up to 3.94 Mm with a resolution of 12.31 km vertically and 20.00 km horizontally. The extension of the 3D box is needed to be able to properly synthesize a spectral-line profile (see Section 5). The box is placed in the upper part of the convection zone with approximately 500 km belonging to the photosphere in 2D and 650 km in 3D, defining the surface layer, h = 0, as the one having a horizontally averaged temperature equal to the solar effective temperature.
The temporal integration is performed using the Strong Stability Preserving Runge-Kutta scheme SSP-RK(3,2) (Kraaijevanger, 1991;Kupka et al., 2012). This is the three-stage scheme of order two with maximum Kraiijevanger radius, which yields the largest possible time step that does not violate the strong-stability condition. Its superior efficiency for the usual accuracies achieved by numerical simulations of stellar granulation has been shown by Grimm-Strele, Kupka, and Muthsam (2015). As with any other explicit time integration method the Courant-Friedrichs-Levy (CFL) criterion is ensured by proper limiting of each time step. For the SSP(3,2) scheme, all CFL conditions for each of the terms of the hydrodynamical equations solved by ANTARES are individually checked and ensured by resizing the time step if needed. The Courant numbers of 0.4 for discretization of advection-and diffusion-related terms used in this check are the default values derived from numerical-stability analysis including a safety factor, but they can be lowered by the user, if needed (see also Kupka et al. (2012) and Grimm-Strele, Kupka, and Muthsam (2015)). Furthermore, we use the diffusion approximation for the lower 70 % of the box instead of angular-ray-integration-based calculations since the plasma is optically thick and the diffusion limit holds. This limit is fixed at a geometrical height that ensures that the temperature ( 12,000 K) is high enough and radiative flux ( 1%) is small enough at each grid point during the relaxation and at all later points in time that the diffusion approximation is always fulfilled at that height. The choice of inclinations for the angular-ray integration is done according to the Gauss-Radau quadrature along the polar angle and with an equidistant node distribution along the azimuthal angle with 18 rays for our 3D simulations and 10 rays for the 2D simulations.
The starting model for our simulations is the standard solar model known as "Model S" (Christensen-Dalsgaard et al., 1996). This is a 1D (i.e. vertical stratification) model of the Sun containing the required physical quantities. It uses the same equation of state (EOS) (Rogers, Swenson, and Iglesias, 1996), opacities Ferguson et al., 2005), and chemical composition (Grevesse and Noels, 1993) as ANTARES (see . For the non-gray radiative transfer, the Kurucz opacity distribution functions (ODFs) (Kurucz, 1992(Kurucz, , 1993 are used for line-opacity computation. ANTARES replicates this starting model over columns to have the initial state for multi-dimensional simulations.

Comparison of Radiative Transfer Solvers
As the main goal of the article is to improve the accuracy of the radiative-transfer solver in ANTARES, in this section we present a comparison of snapshots obtained with the different radiative-transfer solvers. To save computational costs for this comparison we have run the simulations in 2D geometry.
To understand which of our modifications is the most significant one, we have run four types of simulations. The first one uses the original ANTARES radiative-transfer solver with linear integration for optical depth and linear extrapolation of the thermodynamical quantities at the boundaries. It also contains checks for the quantities to be positive. Hereafter, we call this simulation RT1. To understand how strongly the different extrapolation schemes at boundaries affect the radiation field in each simulation, we present the second simulation (RT2) where the same solver with the modified extrapolation at the upper boundary is used. The third simulation (RT3) contains all of our modifications, namely "BézierCP2" radiativetransfer solver, "BézierCP1" scheme for the optical depth calculation, and the modified upper boundary, and RT4 is simulated with linear extrapolation of the upper boundary but using the newly implemented radiative-transfer solver. For all simulations we keep the same initial conditions. We evaluate results for the different cases at the end of respective ANTARES runs that we evolved for a relatively short time (five minutes of solar time) in order for them to have physical conditions not too far away from each other. The flux difference between the different schemes at that resolution in 2D corresponds to a difference in effective temperatures of about 30 K.
The radiative-heating rate is the one quantity for which the accurate radiative-transfer solution is needed, as this is the quantity that appears in the dynamical equation for the temporal evolution of energy solved by ANTARES. Therefore, first we examine Q rad computed in RT1 and RT3 simulations. As is seen in Figure 4, the major difference between the two simulations is seen at the surface of granules, where RT3 shows smoother and numerically less noisy surface structures compared to RT1.
For a more detailed analysis we choose another quantity, radiative flux (see Equation 5). Even though we do not use Equation 7 for the radiative-heating rate computation, the radiative flux is an interesting diagnostic quantity used for visualization of the effects that are hidden in Q rad , such as ray artifacts, but that can also cause some instabilities. In Figure 5 we present the radiative flux calculated with the three radiative-transfer schemes RT1, RT2, and RT3. Similar to Q rad , the major difference between RT1, RT2, and the RT3 simulations is seen at the surface of granules, where RT3 shows smoother and numerically less noisy surface structures compared to RT1 and RT2.
Usually, to save computational costs, all 3D simulations are done with only 18 rays and for all 2D simulations only ten rays were used along which the radiative-transfer equation is solved (see Section 4.1). Because of the limited number of rays, the simulations are relatively fast, but this also introduces some artifacts in the radiative flux in the atmospheric layers (see Figure 5 and 6). To reduce these artifacts either more rays or a more accurate Figure 4 The radiative-heating rate simulated with RT3 (a) and RT1 (b). Note that at a height of ≈ −200 km the artifacts in the downflows forming around a horizontal position of ≈ 4500 km are not present in panel a while being clearly visible in panel b. Also for the two downflows, an improvement is found for the same altitude level. radiative-transfer solver are needed. Figure 6 shows radiative-flux differences in simulations with different radiative-transfer solution, namely using the same modified solver but different top boundary extrapolation, RT4 -RT3 (Figure 6a), and different solvers but the same boundary extrapolation RT4 -RT1 (Figure 6b). It can be seen that the difference between the two solvers is much more prominent than the variation in fluxes due to different boundary extrapolations.
These comparisons clearly reveal the differences between the solvers. However, it is not clear which simulations create more artifacts in the atmospheric layers. In order to clarify this, we enlarge the same region for the simulations done with RT1 and with RT3 (see Figure 7). From a visual comparison of the enlarged figure, it is seen that our new implementation RT3 has fewer artifacts and is more accurate than RT1.

Comparison of Temperature Stratification
The output of ANTARES simulations is a 3D solar or stellar model atmosphere, which is needed in particular to synthesize spectral lines and therefore to interpret fluxes and derive individual characteristics of stars, such as chemical abundances. The 3D simulations with RT1 and RT3 run 55 and 75 minutes of solar time, respectively, during the relaxation phase, and are continued for an additional five minutes of solar time for the production run. We take the latest snapshot for each of the runs and compute one of the products of the radiativetransfer scheme, such as effective temperatures [T eff ] for both simulations, which is 5760 K for RT1 and 5782 K for RT3. The difference of effective temperatures between the two runs for these particular snapshots is around 20 K, which may be due to both the different numerical methods and to taking them at slightly differently relaxed states. In the previous subsection, we showed that our new implementation of the radiative-transfer solver leads to a smoother surface structure and to fewer artifacts in the radiative flux in optically thin layers. The new implementation also affects the temperature stratification, which is one of the important properties of model atmospheres that affect spectral synthesis. In Figure 8 we present an example of the temperature stratification in our 3D simulation cube obtained with RT3 and used for spectral-line synthesis in Section 5. In order for granules to be visible, we show the upper layers in the cube as partially transparent.
We compare the models in 1D geometry because it provides a simple approach while giving sufficient insight. The derivation of a 1D model atmosphere is possible by averaging the 3D model in the horizontal directions at different reference depth-points, e.g. constant heights, column mass, or optical depth. The averaging method should be selected carefully, depending on the application, as different techniques provide quite different models. For more details on the effect of different averaging techniques, we refer the reader to Magic et al. (2013b). As in this article we aim at comparing the temperature structures derived from RT1 and RT3 simulations, we apply a simple geometrical average using the arithmetic mean, which preserves the conservation properties of energy as well as hydrostatic equilibrium.
In Figure 9 we present the comparison of the temperature-height dependence for different solar models and their differences. In order to understand how good the averaged ANTARES models are, in addition to showing the 1D models derived on the basis of the 3D solution for the original RT1 and newly implemented RT3 radiative-transfer solvers, we compare them to other 1D solar models, such as ATLAS (Kurucz, 1992), the FALC semi-  empirical model for the quiet Sun, and Model S, which is also the starting model for the ANTARES simulations. The various stratifications of the models were aligned to each other at the height where T = T eff . In the FALC model, the temperature stratification is derived from observations considering the height dependence of line formation or center-to-limb variations of the continuum. Figure 9 demonstrates that the differences are not relevant for layers where T > 5000 K for almost all models, but quite large differences appear in the upper part of the atmosphere. The ATLAS model deviates from other models in the deep layers, which are hotter than  7000 K where the convection plays a role. It can be due to small mixing-length parameter (α = 1.25) in 1D model atmosphere computations. As none of the models except for FALC include the physical processes needed to form a chromosphere, we cannot trust any of the models above the temperature minimum (500 km). However, the structure below 400 km should be reliable. The reason to extend the 3D simulation cube beyond that range is, as already mentioned by Asplund et al. (2000), to reduce the influence of the upper boundary conditions on the bulk of the simulations. We see that the deviations between RT1 and RT3 are the largest at higher photospheric layers with T < 5000 K; however, to clarify which of the two results is better, we compare them with the other models. We find that RT3 generally agrees better with FALC and ATLAS than RT1.
To understand how good our simulations are in the deeper layers (below the optical surface) it is necessary to compare them with a solar model that has the correct adiabat in the solar interior and not just with atmospheric models. This is ensured by comparing the ANTARES simulation to Model S, because the latter reproduces the depth of the solar convection zone in addition to the solar radius and hence its adiabat agrees with helioseismic observations. Consequently, this match confirms that the layers of the simulation located below the super-adiabatic peak have the correct temperature structure. In Figure 9, we show that both of our models have a very good agreement with Model S in the deeper layers. For additional comparison of ANTARES simulations with other 3D simulations we refer the reader to Kupka (2009).
Another quantity that we want to study in this section, and which plays an important role in spectral-line synthesis, is the temperature gradient. Similar to Figure 9, in Figure 10 we present four panels, one with a direct comparison of considered models and three additional panels with a difference of the gradient between the models. The largest difference between RT1 and RT3 is in the near-surface layers where the continuum forms and where the gradient is largest. However, to understand which model yields the more consistent temperature gradient with other models we show the difference with the gradient derived from ATLAS, FALC, and Model S. It is clearly seen that in these layers the deviation of ATLAS -RT3, FALC -RT3, and Model S -RT3 is less than the difference with RT1. This result allows us to conclude that improving the radiative-transfer solver in hydrodynamic simulations leads to an improvement of the temperature structure and its gradient.

Application: Spectral Line Profiles
In the previous section, we have shown that the implementation of the more accurate radiative-transfer solver improves the appearance of typical structures in our 3D hydrody- Table 1 Atomic parameters of the chosen spectral line applied to the line synthesis are taken from VALD-3 Database (Piskunov et al., 1995;Kupka et al., 1999;Ryabchikova et al., 2015).

Wavelength [Å]
Excitation namical convection simulations. As an important astrophysical application in this section, we describe how these improvements affect the profile of photospheric spectral lines. As in this article we do not aim at deriving any astrophysical parameters from the considered spectral lines but mainly focus on visual comparison of RT1 and RT3 with observations, we limit ourselves to only one snapshot and do not discuss a temporal average of many snapshots.

Spectral-Line Selection
Given the magnetic-field-free property of our simulations, we select non-magnetic or weakly magnetic Fe I spectral lines in the visible spectral range, with a laboratory wavelength of λ = 5569.62 Å, 5577.02 Å, and 5584.76 Å for consistency. We expect spectral synthesis with these lines to be reliable. We provide the adopted atomic parameters of these lines in Table 1. The spectral line λ = 5584.76 Å is compared with different computations in the next subsection. To understand the selection of the spectral line, being the least magnetically sensitive of the three we selected, it is worth showing where it forms. A useful diagnostic tool for studying the spectral-line formation heights is the line-contribution function. It displays the contribution to the intensity at different heights. We highlight the formation heights of the selected spectral line through a visual representation of its contribution functions at the disk center, as shown in Figure 11. It is seen that the chosen line forms in layers of the solar atmosphere that are within the upper vertical extension of our chosen ANTARES simulation box and are between about −150 km (thus, a bit below the optical surface) and up to about 400 km above it. This is the region where both the temperature stratification and the temperature gradient are more reliable for RT3 than for RT1 and therefore we expect to obtain a better agreement of the spectral line synthesized from the RT3 model atmosphere with observations.

Spectral-Line Synthesis
We compare the synthetic spectra at the disk center from the ANTARES simulation from two runs, one using RT1 and the second using RT3 as radiative-transfer solver in ANTARES, to understand how the accuracy of the radiative-transfer solver in hydrodynamic simulations affects the spectral-line synthesis. The model atmosphere for the synthesis consists of the upper layers of one 3D ANTARES simulation snapshot. To compute the high-resolution spectral line we use the Merged Parallelized Simplified ATLAS code (MPS-ATLAS: Witzke et al., 2021) developed from the original ATLAS code (Kurucz, 1970). For consistency, we use the same abundances (Grevesse and Sauval, 1998) for all spectral-line synthesis calculations as in the OPAL tables we employed in our ANTARES simulations (see Section 4.1). As MPS-ATLAS allows the use of different line lists, we chose the newest Vienna Atomic Line Database (VALD-3) line lists (Piskunov et al., 1995;Kupka et al., 1999;Ryabchikova Figure 11 Contribution function of the λ = 5584.76 Å spectral line in the ANTARES RT3 horizontally averaged model atmosphere. et al., 2015) with continuum co-added from ATLAS to synthesis high-resolution spectral line. Note that to take into account line blanketing in ANTARES simulations, the line opacities from the ATLAS opacity distribution function (ODF) package (Castelli, 2005) including continuum opacities of ATLAS sorted into four opacity bins are used (for more details, see Muthsam et al., 2010). The spectral-line broadening due to velocity fields was introduced into the MPS-ATLAS calculations by considering a constant microturbulent velocity of 0.85 km s −1 and a macroturbulence of 1.5 km s −1 and convolving the averaged spectralline profile with the corresponding Gaussian profile (see details in Asplund et al., 2000). As we only aim at a visual comparison of our simulations with disk-center observations, we synthesize the 3D spectral lines for the center of the solar disk. Note that in MPS-ATLAS a different radiative-transfer scheme, modified Feautrier (Lester and Neilson, 2008), is used.
As MPS-ATLAS takes into account only constant values of micro-and macroturbulent velocities, we additionally compute the lines from the 3D simulation cube using the SPECTR-3D post-processing tool (Piskunov, 2021), where the full velocity field from the ANTARES snapshot is taken into account in the radiative transfer. SPECTR-3D uses the same equation of state as in ANTARES and the same VALD-3 line lists and continuum that we use in MPS-ATLAS. In order to study how the velocity field affects the spectral-line synthesis, we compute the spectrum only for the 3D simulation run with RT3 and compare with the MPS-ATLAS synthesis for the same ANTARES run. Note that the radiative-transfer solver used in this tool is the Bézier radiative-transfer scheme (de la Cruz Rodríguez and Piskunov, 2013). For a direct comparison with MPS-ATLAS hence some additional uncertainties could result from using different numerical methods, although at the chosen resolution they should be acceptably small (see the latter reference regarding tests of the influence of numerical resolution, in particular their Figures 2 and 3).
To understand how well the synthetic spectral-line profiles reproduce the observations, we compare them with the solar disk-center intensity atlas observed with the Fourier-Transform Spectrometer (FTS) at the McMath-Pierce telescope (Neckel, 1999), also dubbed the "Hamburg atlas". We prefer the FTS solar atlas to the "Photometric Atlas of The Solar Spectrum from λ3000 to λ10000" published by the University of Liège (Delbouille, Roland, and Neven, 1973) because of the higher spectral resolution of λ/ λ ≈ 500,000. As was shown by Doerr, Vitas, and Fabbian (2016), where the authors carry out a detailed comparison of FTS and Liège atlases, the spectral resolution of the Liège atlas is between two and Figure 12 Comparison of the λ = 5584.76 Å spectral line from MPS-ATLAS and SPECTR-3D synthesis and in the Hamburg solar atlas at the solar disk center (Neckel, 1999). six times lower. The wavelength shift due to the Doppler effect from the relative Sun-Earth movement is already corrected in the observational data. We additionally removed the redshift (633 m s −1 ) of the observed spectral-line profiles to take into account the relativistic effect caused by the solar gravitational field.
In Figure 12 we present the comparison of our synthetic spectral lines with observations. First of all, we want to point out that the MPS-ATLAS line computed for the RT3 run reproduces the observations slightly better in the line core than the one for the RT1 simulation, showing that the new implementation improves the model atmosphere temperature stratification. The MPS-ATLAS spectral-line synthesis has a blend on the blue side and therefore it does not match the observations. This blend is absent in the SPECTR-3D spectral-line synthesis. On the red side the bad match is due to the symmetric nature of the micro-and macroturbulence approximations of the MPS-ATLAS calculations.
The second important point that we want to stress here is that the SPECTR-3D spectral line matches the observations better. The center of the line is too weak compared to the observations, a discrepancy that could be improved by increasing the height of our simulation box. However, the width of the line nicely matches the observed one. This implies that the velocity field from our ANTARES simulations is correct. In Figure 13 we present a larger spectral range of disk-center intensity observations and SPECTR-3D computations with additional panels zooming in on a strong and a weak spectral line. It demonstrates that the spectrum computed from 3D ANTARES simulations that are based on the RT3 scheme reproduces the observed spectrum. The weaker line that we show as an example matches observations while the strong line to which we zoom in is too broad. We want to note here that we do not tune any abundances to have better agreement with observations. This line is slightly sensitive to the magnetic field (see Table 1) and most likely there is a blend in the red wing in the SPECTR-3D calculation, that seems absent in the observations, pointing to an atomic-data problem. Additionally, this line with equivalent width of ≈ 130 -140 mÅ is too strong for our simulation box as its core forms at the top of the box. Such a strong line can reliably be modelled only by 3D RHD simulations that properly account for higher-lying layers. It would thus require a proper modeling of the chromosphere, which is not within the scope of this study. For example, Asplund et al. (2000) have studied in detail many iron lines, selected to all have an equivalent width of no more than 100 mÅ. The effects just mentioned all influence the resulting shape of the spectral-line profile and will be studied in a separate article. As was already mentioned, the predicted profiles of stronger lines Figure 13 Comparison of SPECTR-3D synthetic spectral lines (red) and observational data from the Hamburg solar atlas (Neckel, 1999) (blue) at the disk center.
would benefit from using a simulation box extending higher into the upper photosphere and chromosphere.

Discussion
Accuracy and stability of radiative-transfer solvers are crucial ingredients of 2D and 3D HD simulations used in astrophysics. Typically, the grid resolution of such simulations is low in the sense that there are regions in the flow where the optical depth, even in the sense of the Rosseland average, changes rapidly within a few grid cells in a depth range of order unity, i.e. where the emergent intensity arises. Thus, stable and accurate integrations of the optical depth and the source function are needed to ensure reliable calculations of the radiative cooling and heating rates. The Bézier-spline-based method by de la Cruz Rodríguez and Piskunov (2013) provides a solution to this challenging problem from a numerical point of view. However, in the context of 2D and especially 3D HD, scaling properties of the simulation code with respect to its performance on large-scale parallel systems are important and for domain-based parallelization concepts, such as the one used in ANTARES, locality of the data required for the calculations is simplifying efficient parallelization.
Indeed, when implementing the full method of de la Cruz Rodríguez and Piskunov (2013) into a 3D hydrodynamic simulation program, a problem arises. In order to implement the full Bézier technique a fourth point (k − 2) would have to be added. On the other hand, its simplification to a "BézierCP1" scheme decreases the accuracy of the Bézier method, whereas the full method by de la Cruz Rodríguez and Piskunov (2013) is expected to lead to better stability and accuracy. However, the full method necessarily decreases the speed of calculations as an additional interpolation for the fourth point between nodes has to be performed. This is also a less local procedure with respect to the data the algorithm needs access to. For this reason we have suggested here the simplifications provided by the "BézierCP2" scheme and "BézierCP1" scheme for the integration of the source function and the optical depth, respectively. This way, accuracy is gained in comparison with purely linear (trapezoidal rule) based integration and also the stability is enhanced when comparing the new method to the parabolic one originally implemented into the ANTARES code.
In comparing synthetic solar spectra computed on the basis of our new 3D HD simulations with observations, we have also provided another, successful demonstration of the physical realism of ANTARES-based models of the solar photosphere. Regarding further possible applications of synthetic spectra, for example to solar-and stellar-irradiance and -variability studies, we point out that one should remain aware of the challenges and uncertainties involved, such as spectrum calibration, opacity selection, and so on; see for example Criscuoli et al. (2020). In particular, this is true in the context of solar-irradiance reconstructions, given the uncertainties in the needed input data for example opacities and atomicdata sources. In this same topical collection, Criscuoli et al. (2020) highlight that different radiative-transfer codes and methodologies produce different results on synthetic radiative fluxes of quiet and magnetic solar-surface regions, especially concerning contrast of magnetic features. Matching very accurately across a broad range of wavelength regions the flux spectrum observed to emerge from quiet and magnetic features on the solar atmosphere remains a significant challenge.

Summary
In this article, we have presented the new radiative-transfer solver implemented in ANTARES 3D radiation hydrodynamical simulations. The 3D simulations are applied to a range of tasks in solar and stellar physics. Most of the existing 3D HD/MHD codes, such as STAGGER, MURaM, or CO 5 BOLD, implement the radiative-transfer solver in a very simple way (linear interpolation of the source function and optical depth) to reduce the computational costs. The linear integration is stable; however, at lower resolutions it cannot reproduce the smooth structure of granules. We have shown that the new solver predicts smoother surface structures for the simulated granules and decreases the number of artifacts in the photosphere and optically thin part of the simulation cube, where radiative transfer is the main mechanism of energy transport. Moreover, because the Bézier interpolation is monotonic and additional flux limiters are no longer necessary, the speed of the computations with the new radiative-transfer solver is almost the same despite the higher order of the scheme, and hence it does not slow down the simulations significantly. Therefore, from now on ANTARES radiation hydrodynamical simulations will be done with the RT3 method described in this article.
We studied how the averaged 1D temperature structure and its gradient (both derived from 3D simulations with the original and the newly implemented radiative-transfer solver) are stratified. We compared them with the 1D ATLAS solar model atmosphere and the FALC semi-empirical model atmosphere for the quiet Sun, and we showed that the horizontally averaged temperature of the simulation carried out with RT3 shows better agreement with the ATLAS and FALC temperature structure. The atmospheric temperature gradient is closer to the gradient in ATLAS, FALC, and Model S, which is the starting model of ANTARES simulations used here. Therefore we conclude that a more accurate radiative-transfer solver in radiation hydrodynamical simulations leads to a more accurate stratification of the temperature and of its gradient at typical resolutions of HD simulations of solar and stellar surface convection.
We present one of the applications of 3D hydrodynamical simulation, namely spectralline synthesis. We synthesized magnetically insensitive Fe I spectral lines using SPECTR-3D and MPS-ATLAS post-processing tools applied to the 3D ANTARES model atmosphere obtained with the original and the new radiative-transfer solvers. To understand which approach performs better we compared the synthetic spectra with the Hamburg FTS atlas of the solar spectrum. We show that our new implementation improves the fit of spectral-line profiles, which could be useful for solar abundances derivation. The weaker lines, which form in the lower atmosphere, are very well reproduced. Somewhat stronger spectral lines could be modeled properly when increasing the vertical extent of the simulation cubes. A comparison of the profiles of the spectral lines synthesized from our new implementation, RT3, with the one from the original radiative-transfer solver, RT1, shows that the more accurate solver in HD simulations displays better agreement with observations.