Issue 
A&A
Volume 675, July 2023



Article Number  A160  
Number of page(s)  26  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202346341  
Published online  18 July 2023 
Opacity for realistic 3D MHD simulations of cool stellar atmospheres
^{1}
Instituto de Astrofísica de Canarias,
38200
La Laguna, Tenerife,
Spain
email: andrperdomo@gmail.com
^{2}
Departamento de Astrofísica de la Universidad de La Laguna,
38200
La Laguna, Tenerife,
Spain
^{3}
Isaac Newton Group of Telescopes,
Apartado de Correos 321,
38700
Santa Cruz de La Palma, Canary Islands,
Spain
^{4}
Steward Observatory, University of Arizona,
Tucson,
USA
Received:
7
March
2023
Accepted:
5
June
2023
Context. Realistic threedimensional timedependent simulations of stellar nearsurface convection employ the opacity binning method for the efficient and accurate computation of the radiative energy exchange. The method provides several orders of magnitude of speedup, but its implementation includes a number of free parameters.
Aims. Our aim is to evaluate the accuracy of the opacity binning method as a function of the choice of these free parameters.
Methods. The monochromatic opacities computed with the SYNSPEC code were used to construct opacity distribution function (ODF) that was then verified through detailed comparison with the results of the ATLAS code. The opacity binning method was implemented with the SYNSPEC opacities for four representative cool mainsequence stellar spectral types (F3V, G2V, K0V, and M2V).
Results. The ODFs from SYNSPEC and ATLAS show consistent results for the opacity and bolometric radiative energy exchange rate Q in the case of the F, G, and Ktype stars. Significant differences, coming mainly from the molecular line lists, are found for the Mtype star. It is possible to optimise a small number of bins to reduce the deviation of the results coming from the opacity grouping with respect to the ODF for the F, G, and Ktype stars. In the case of the Mtype star, the inclusion of splitting in wavelength is needed in the grouping to get similar results, with a subsequent increase in computing time. In the limit of a large number of bins, the deviation for all the binning configurations tested saturates and the results do not converge to the ODF solution. Due to this saturation, the Q rate cannot be improved by increasing the number of bins to more than about 20 bins. The more effective strategy is to select the optimal location of fewer bins.
Key words: opacity / radiative transfer / equation of state / stars: atmospheres
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Radiative transfer is one of the key factors governing the structure and dynamics of stellar atmospheres and it has to be taken into account with special care when constructing a model atmosphere (Hubeny & Mihalas 2014). Detailed interactions between radiation and matter are quantitatively described by the opacity, that is by the absorption coefficients of all contributing (boundbound, boundfree, and freefree) processes due to all relevant particle species present in the atmosphere (Rutten 2003). The complexity of the opacity function is enormous as it is visible from the plethora of atomic and molecular spectral lines in the spectra of latetype stars (Kurucz 2002, 2017) and its evaluation is computationally expensive. In combination with the already computationally demanding solution of the radiative transfer equation (RTE), the problem of calculating detailed monochromatic radiative losses and gains in the atmosphere, accounting explicitly for ≈10^{9} lines and ≈10^{6} wavelength points, is difficult even in onedimensional (1D) modelling, and virtually impossible in threedimensional (3D) timedependent numerical simulations. Two important tricks are available, however, that reduce the computational effort by many orders of magnitude.
First, under the assumption of local thermodynamic equilibrium (LTE), opacity is a function of chemical composition, the turbulent velocity, wavelength, and two independent thermodynamic variables (for example, temperature T and mass density ρ); it can therefore be precomputed and stored into lookup tables. Even in the situations where some atoms are treated out of LTE, the opacity of the LTE background can be precomputed without the contribution of the nonLTE (NLTE) lines, which are computed explicitly. The opacities are computed as byproducts in the codes and projects for modelling stellar atmospheres: the ATLAS family of codes (Kurucz 1970; Castelli & Kurucz 2003; Sbordone et al. 2004; Witzke et al. 2021), MARCS (Gustafsson et al. 1975, 2008), PHOENIX (Gustafsson et al. 1975, 2008), SYNSPEC/TLUSTY (Hubeny & Lanz 2011; Hubeny et al. 2021), and others.
Second, in the energy conservation equation of the system of magnetohydrodynamic equations, the radiative energy exchange is accounted for by a wavelengthintegrated source term Q. The distribution of radiative heating and cooling with wavelength is irrelevant, at least as long as the onefluid magnetohydrodynamics (MHD) approach is used. Therefore, instead of solving monochromatic RTE and then integrating the results, an approximate solution may be obtained by first combining or integrating the opacities in wavelength ranges and then solving the RTE for a number of representative opacity values. This idea is implemented in the opacity distribution function (ODF) defined by Strom & Kurucz (1966). In their approach the monochromatic opacities are divided into hundreds of wavelength segments. The opacities in each segment are then sorted by magnitude to obtain a continuous monotonie distribution that is then discretised by averaging over several (around 10) substeps. In that way the RTE can be solved for several thousands of ODF values instead of for millions of monochromatic opacities.
An important underlying assumption in the ODF approach is that there are no substantial wavelength shifts in the spectrum along the different heights in the atmosphere. There are situations in which this assumption fails (see chapter 17.6 in Hubeny & Mihalas 2014), especially for plasmas with velocity fields. Nevertheless, even in moving media, if we are not interested in detailed spectra and if the velocities are not larger than the thermal velocity, we can still solve the RTE using ODFs.
Precomputed ODF tables produced by the ATLAS code are widely used in stellar modelling and for computing emergent spectra from the models (see Kurucz 2014). The ODF method was recently optimised and generalised by Cernetic et al. (2019). Nevertheless, not even the speedup enabled by ODF is sufficient when atmospheres are modelled as timedependent 3D systems. A solution, known in the literature as the opacity binning (OB) method, was initially proposed by Nordlund (1982). In this method the opacities are grouped according to the optical depth in a representative atmosphere at which the monochromatic optical depth reaches unity. Very few groups (or bins) are sufficient to approximate fairly accurately the correct solution for the radiative energy exchange rate, although the explicit wavelength dependence of the opacity is eliminated. The number of groups and the distribution of the optical depth separators between them are free parameters of the method, while four groups are typically used, roughly representing the continuum, weak, intermediate, and strong lines. The OB method opened the door to the efficient implementation of realistic 3D simulations that have revolutionised our understanding of the physics of stellar atmospheres. The method was further developed by Ludwig (1992), Ludwig et al. (1994), Vögler (2003), and Trampedach et al. (2013). Skartlien (2000) developed a variant that includes the effects of scattering.
The wavelength dependence of the opacities is disregarded in the OB method. Vögler et al. (2004) and Ludwig & Steffen (2013) noted that this approximation does not converge to the correct solution in the limit of an infinite number of bins. Several other groups experimented with variations to the OB method. Caffau et al. (2008) used 12 bins in their study of the oxygen abundance, while Pereira et al. (2013) and Magic et al. (2013) introduced a variant of the method whereby the opacities are sorted by wavelength into several groups before sorting by optical depths. However, there is no unique criterion for choosing the number and distribution of the bins. The distribution of the opacities varies significantly with the effective temperature and chemical composition of a star and it may thus be expected that different distributions of bins are optimal for different stars. Our aim in this paper is to test the OB method for four cool mainsequence stars (F3V, G2V, K0V, and M2V) and to design a strategy for finding optimal setups. Our analysis is based on monochromatic opacities computed with the SYNSPEC code. The code (and its modelling counterpart TLUSTY) was recently significantly upgraded and equipped with an uptodate list of atomic and molecular spectral lines necessary for the modelling of cool stellar types (Hubeny et al. 2021). The code is publicly available as open source^{1} and supplemented with a Python wrapper (Allende Prieto et al. 2022). As we intend to use the code to prepare customised opacities for 3D simulations of nearsurface convection in cool stars (Perdomo et al., in prep.) with the MANCHA code (Khomenko et al. 2017, 2018), flexibility in selecting opacity contributors is of particular importance for us.
In Sect. 2 we briefly introduce a set of the ID models that represent our stars and specify the radiative transfer details. Computation of the monochromatic opacities using the SYNSPEC code are described in Sect. 3. In Sect. 4 the essentials of the ODF method are summarised: ODFs are constructed from the monochromatic opacities computed with SYNSPEC, and the results are compared with the results of the ATLAS code. The same data set is used to study strategies for OB in Sect. 5: we first test the importance of the location of the bin separators in optical depth (Sect. 5.1), then we explore the influence of the number of bins in optical depth (Sect. 5.2), and, finally, we study two different strategies for combined binning in optical depth and wavelength (Sect. 5.3). Our conclusions are presented in Sect. 6.
Fig. 1 Temperature stratification of the F3V (green), G2V (yellow). K0V (orange), and M2V (dark red) 1D hydrostatic models used in this work. The same colour coding is used throughout the paper to label results for the four stellar types. The grid used to solve the RTE is marked by dots. 
2 Method
2.1 Model atmospheres
To test the OB for a series of main sequence spectral types, a set of 1D stellar atmospheric models is used. These atmospheres are computed assuming hydrostatic equilibrium and using the mixing length theory to account for the convective energy transfer. The models are initiated with radiation in the grey approximation and corrected to have representative temperature gradient at the top of the atmosphere by scaling the Hopf function q(τ) (see chapter 17 from Hubeny & Mihalas 2014, and references 507– 510 there) of the HarvardSmithsonian Reference Atmosphere (Gingerich et al. 1971) with the effective temperature of each star, following the expression ${T}^{4}\left(\tau \right)=\frac{3}{4}{T}_{\text{eff}}^{4}\times \left[\tau +q\left(\tau \right)\right]$. The models are shown in Fig. 1 and are computed for spectral types F3V, G2V, K0V, and M2V, with effective temperature (K) and logarithm of gravity (cm s^{−2}) of [6890, 4.3], [5780, 4.4], [4855, 4.6], [3690, 4.8], respectively. In the computation of the models we used the Rosseland opacities consistent with the opacities used in the present paper (see Sect. 3). More details about the models will be published elsewhere.
2.2 Radiative transfer computation
The RTE in a planeparallel static atmosphere in LTE without considering line scattering is
$$\mu \frac{d{I}_{v}\left(z,v\right)}{dz}=\rho \left(z\right){\varkappa}_{v}\left(z\right)\left[{B}_{v}\left(z\right){I}_{v}\left(z\right)\right],$$(1)
where I_{v}(z, ν) is the specific intensity (erg s^{−1} cm^{−2} Hz^{−1} ster^{−1} ) at frequency ν in the direction µ and at the geometrical height z in the atmosphere; B_{v} is the Planck function; ρ(z) is the mass density; and ϰ_{ν} (z) is the monochromatic absorption coefficient per unit mass (cm^{2} g^{−1} ; referred to as opacity throughout the text). The opacity is computed as the sum of the continuum opacity (mainly contributions from the freefree and boundfree processes) and the line opacity (boundbound individual transitions included in the line lists for atoms and molecules),
$${\varkappa}_{v}={\varkappa}_{v}^{\text{cont}}+{\varkappa}_{v}^{\text{lines}}.$$(2)
We restrict ourselves to the LTE case that is valid only in deep layers where J ≈ B. Therefore our treatment of the coherent scattering terms (see Appendix A) in the total opacity is only approximate (e.g. Hayek 2010).
We solve the RTE in the vertical direction (µ = ±1). At each frequency ν, the formal solution for outward and inward intensity (I^{+}, I^{−}) in LTE is:
$$\begin{array}{lll}{I}_{v}^{\pm}\left({z}_{i}\right)\hfill & =\hfill & {I}_{v}^{\pm}\left({z}_{i\mp 1}\right)\mathrm{exp}\left(\text{\Delta}{\tau}_{v}^{i\mp 1}\right)+{\displaystyle {\int}_{0}^{\text{\Delta}{\tau}_{v}^{i\mp 1}}{B}_{v}\left(t\right)\mathrm{exp}\left(t\text{\Delta}{\tau}_{v}^{i\mp 1}\right)\text{d}t}\hfill \\ \hfill & \equiv \hfill & {I}_{v}^{\pm}\left({z}_{i\mp 1}\right)\mathrm{exp}\left(\text{\Delta}{\tau}_{v}^{i\mp 1}\right)+\text{\Delta}{I}_{v}^{\pm}\left({z}_{i}\right),\hfill \end{array}$$(3)
where i is the height index for each point in the atmosphere, being zero at the bottom and increasing upwards; ${\tau}_{v}^{i}$ is the optical depth at the height with index i (dτ_{v}(z) = ρ(z)ϰ_{ν} (z)dz); and $\text{\Delta}{\tau}_{v}^{i+1}={\tau}_{v}^{i+1}{\tau}_{v}^{i}$ and $\text{\Delta}{\tau}_{v}^{i1}={\tau}_{v}^{i}{\tau}_{v}^{i1}$. To solve the RTE, we apply the short characteristics method (Olson & Kunasz 1987), using the linear approximation of the Planck function:
$$\text{\Delta}{I}_{v}^{\pm}\left({z}_{i}\right)={\psi}_{0}{B}_{v}\left({z}_{i}\right)+{\psi}_{1}{B}_{v}\left({z}_{i\mp 1}\right).$$(4)
For each frequency ν we compute the coefficients for the local and the upwind points
$$\begin{array}{l}{\psi}_{0}=1\frac{1}{\text{\Delta}{\tau}_{v}}+\frac{1}{\text{\Delta}{\tau}_{v}}\mathrm{exp}\left(\text{\Delta}{\tau}_{v}\right),\hfill \\ {\psi}_{1}=\frac{1}{\text{\Delta}{\tau}_{v}}\frac{\text{\Delta}{\tau}_{v}+1}{\text{\Delta}{\tau}_{v}}\mathrm{exp}\left(\text{\Delta}{\tau}_{v}\right).\hfill \end{array}$$(5)
Since the method is formulated in a way that it is symmetrical along a ray direction, the coefficients are the same for the upward and inward intensities, and the only change needed is the specific Δτ_{ν} in each case ($\text{\Delta}{\tau}_{v}=\text{\Delta}{\tau}_{v}^{i1}$ for the upward intensities; $\text{\Delta}{\tau}_{v}=\text{\Delta}{\tau}_{v}^{i+1}$ for the inward ones). This first order shortcharacteristics scheme would not be sufficiently accurate if scattering was taken into account (i.e. if the formal solver was used in an iterative solution).
The mean intensity J_{v} (erg s^{−1} cm^{−2} Hz^{−1} ster^{−1}) and the flux ℱ_{v} (erg s^{−1} cm^{−2} Hz^{−1}) can then be calculated for a vertical ray from the known I_{v}:
$${J}_{v}\left({z}_{i}\right)=\frac{1}{2}{\displaystyle {\int}_{1}^{1}{I}_{v}\left({z}_{i}\right)\text{d}\mu \equiv \frac{1}{2}\left({I}_{v}^{+}\left({z}_{i}\right)+{I}_{v}^{}\left({z}_{i}\right)\right),}$$(6)
$${\mathcal{F}}_{v}\left({z}_{i}\right)=2\pi {\displaystyle {\int}_{1}^{1}{I}_{v}}\left({z}_{i}\right)\mu \text{d}\mu \equiv 2\pi \left({I}_{v}^{+}\left({z}_{i}\right){I}_{v}^{}\left({z}_{i}\right)\right).$$(7)
Finally, the monochromatic radiative energy exchange rate Q_{v} (erg cm^{−3} s^{−1} Hz^{−1}), which accounts for the radiation sources and sinks in the energy conservation equation, can be computed either as the divergence of the radiative flux,
$${Q}_{v}^{F}\left({z}_{i}\right)=\nabla {\mathcal{F}}_{v}\equiv \left({\mathcal{F}}_{v}\left({z}_{i+1}\right){\mathcal{F}}_{v}\left({z}_{i1}\right)\right)/\left({z}_{i+1}{z}_{i1}\right),$$(8)
or directly from the mean intensity,
$${Q}_{v}^{J}\left({z}_{i}\right)=4\pi {\varkappa}_{v}\left({z}_{i}\right)\rho \left({z}_{i}\right)\left({J}_{v}{B}_{v}\right).$$(9)
The monochromatic rate is calculated as a weighted combination of the two, using Eq. (4.13) from Vögler (2003; following the suggestion from Bruls et al. 1999):
$${Q}_{v}\left({z}_{i}\right)=\mathrm{exp}\left({\tau}_{v}^{i}/{\tau}_{\text{h}}\right){Q}_{v}^{J}\left({z}_{i}\right)+\left[1\mathrm{exp}\left({\tau}_{v}^{i}/{\tau}_{\text{h}}\right)\right]{Q}_{v}^{F}\left({z}_{i}\right),$$(10)
where τ_{h} = 0.1. As it is explained in Vögler (2003), this smooth transition between ${Q}_{v}^{J}$ and ${Q}_{v}^{F}$ avoids the accuracy problem that ${Q}_{v}^{J}$ has in the optically thick regime (where J_{v} approaches B_{v}) and the errors coming from small variations of the orientation of the flux in the optically thin part of the atmosphere (where flux should be nearly constant), that are enhanced by the gradient in ${Q}_{v}^{F}$.
3 Monochromatic opacity
For the computation of the monochromatic opacity ϰ_{v}, we use synple , a Python wrapper for SYNSPEC . SYNSPEC is a general spectrum synthesis code that has been used to solve the RTE in different astrophysical scenarios, including cool stars (see, for example, Allende Prieto et al. 2003b,a and Palacios et al. 2010; Osorio et al. 2020). The main purpose of SYNSPEC is to synthesise spectra for a given model atmosphere, but it can also be used to compute lookup tables of monochromatic opacities for a given grid of thermodynamic quantities. The code solves the equation of state (EOS), and uses the number densities to compute the continuum and line opacities. The EOS from SYNSPEC is fully consistent with the results of our own EOS solver (Vitas et al., in prep.), which is used to produce the EOS tables for the 3D simulations of stellar atmospheres with the MANCHA3D code (Felipe et al. 2010; Khomenko et al. 2017, 2018; Modestov et al., in prep.).
The EOS in SYNSPEC is solved for a specified grid of temperature T and mass density ρ. The number of nuclei per species (fixed by the abundances) and charge are conserved. The code solves the chemical equilibrium equations to determine the populations for 38 atomic species. For temperatures lower than a certain value (in this work chosen to be 8000 K), neutral atoms and singly ionised atoms are included in the EOS (not taking into account higher ionisations) and the code also considers molecular formation (including 503 molecular species). For higher temperatures, molecules are not included in the EOS, and higher ionisations are computed for the atoms. For more details on the solution of the EOS in SYNSPEC (and other details and the general use of the code), we refer the reader to Section 2.6 from Hubeny et al. (2021), and Hubeny & Lanz (2Ol7a,b).
Once the EOS is solved, the opacity is computed in the same (T,ρ) grid, using the atomic and molecular data from several sources (see Appendix A). The data used by the code can be customised by the user by changing, for example, the line lists for the boundbound transitions. The control over opacity contributors is important for the cooler stars, where omission of certain molecules (e.g. TiO and VO) or incomplete line lists may strongly affect both detailed spectral synthesis and the bulk opacities used for modelling.
For the ingredients listed in Appendix A, we computed a lookup table using synple, with a microturbulence velocity of 2 km s^{−1} and for the solar abundances from Anders & Grevesse (1989). These solar abundances are outdated but preferred in our study since they were used in the computation of the reference ODFs to which they are subsequently compared. The Taxis of the grid has 60 values in the range Τ ϵ [1995, 125 900] K, with constant step Δ log T ≃ 0.03 Κ (in the present text, log ≡ log _{10}). The ρaxis has 35 values in the range ρ ∈ [10^{−13}, 3.2 × 10^{−3}] g cm^{−3} with Δ log ρ ≃ 0.3. The table is computed for wavelengths in the range λ ∈ [20, 9.5 × 10^{4}] nm, with Δ log λ = 10^{−6} (λ in nm). In total there is around 7.7 × 10^{9} data points in the table, stored in an HDF5 file of 230 GB. Another table is computed for the same wavelength grid but with finer temperature and density axes (same number of points, but reduced ranges, Τ ∈ [1995, 19 900] Κ and ρ ∈ [10^{−9},3.2 × 10^{−3}] g cm^{−3}). This table is more suitable for the M2V star, thus reducing the interpolation errors in the RTE solver.
The total opacity (continuum and lines) in the wavelength range from 100 to 2500 nm is shown in Fig. 2. It is interpolated for the temperatures and mass densities of the four 1D stellar models. The opacity clearly changes with different heights in the stellar models and with spectral type. To check where the radiation at each wavelength contributes to the emergent radiation, we mark the geometrical heights at which τ_{ν} = 1 (thick line). This curve includes the contribution from both the line and the continuum opacity. Since the wavelength range is wide, we divided it into subintervals of 2.4 nm and computed the average height in each of them. Similarly, we overplot the formation heights of the continuum (thin line), at which ${\tau}_{v}^{\text{cont}}=1$. The ionisation edges and the typical shape of the boundbound and boundfree H^{−} opacity coefficient are visible for the ${\tau}_{v}^{\text{cont}}=1$ curve. In the UV the line forest elevates the height of the τ_{λ} = 1 curve. The strong spectral lines, especially the hydrogen series, are clearly visible as well as their relative decrease in strength with decreasing effective temperature. In the M star, the τ_{λ} = 1 curve is shaped by millions of molecular transitions.
Following the method from Sect. 2, we solve the RTE to obtain Q_{v}. Once Q_{v} is computed, we can integrate it in frequency using the trapezoidal rule:
$$Q={\displaystyle {\int}_{v}{Q}_{v}\text{d}v\equiv \frac{1}{2}{\displaystyle \sum _{i}\left({Q}_{{v}_{i}}+{Q}_{{v}_{i+1}}\right)\text{\Delta}{v}_{i}}.}$$(11)
The computed Q_{v} is shown in Fig. 3, for the wavelength range and set of stellar models as in Fig. 2. The radiative energy exchange happens mostly around the stellar surface. For decreasing effective temperature, the absolute value of Q_{v} is large for a wider wavelength range and the significance of the contribution to Q_{v} shifts from the visible for the F star towards the IR for the M star. Figure 3 also shows how the relative magnitude and area covered by significant Q_{v} for the heating in the atmosphere (Q_{v} > 0) and the cooling (Q_{v} < 0) changes for the different spectral types. For the F star the maximum heating is around an order of magnitude smaller than the magnitude of the most intense cooling, and peaks in the nearUV and blue wavelengths. For the M star, the heating is comparable to the cooling in magnitude, and ranges from the visible to the nearIR. As in the opacity plot (Fig. 2), the individual contributions to the line opacity from strong transitions are obvious for the hotter stars in the sample. In contrast, they almost disappear for the M star, while small structures due to the presence of the molecules become more evident.
Figure 3 provides a qualitative overview of the monochromatic radiative energy exchange. One must keep in mind that the computations presented here are done in the LTE approximation which is fairly good for most of the photospheric lines, but fails badly for lines that contribute to radiative losses in the chromosphere. However, the method of precomputed opacities is feasible only in the LTE approximation when there is unique mapping between the pair thermodynamic quantities and the computed total opacity. Our primary focus is therefore on the stellar surface (i.e. the top of the convection zone and the lower photosphere), from where the bulk of the radiation escapes the star. This is also the reason why we present the results for the bolometric Q (erg cm^{−3} s^{−1}) and not the ratio Q/ρ (erg g^{−1} s^{−1}), which would put more emphasis on the radiation exchange at lower densities.
Fig. 2 Monochromatic opacity (values shown in the colour–bars) for the F3V (first row), G2V (second row), KOV (third row), and M2V stars (fourth row) versus wavelength (horizontal axis) and geometrical height (vertical axis). Horizontal grey lines: Rosseland optical depths. The zero geometrical height is taken where logτ_{Ros} = 0. Other lines are shifted by Alogτ_{Ros} = 1. Black dashed thin line: height at which the continuum optical depth at each wavelength is unity (${\tau}_{\lambda}^{\text{cont}}=1$). Black solid thick line: heights at which the total optical depth τ_{λ} is unity. Values are averaged over 25 nm intervals. Tick marks on the top indicate wavelengths of the hydrogen series. Hydrogen lines gradually vanish with decreasing effective temperature. The colourbar values of each star range from min log ϰ to max log ϰ. 
Fig. 3 Monochromatic radiative energy exchange rate Q_{v} (values shown in the colourbars) computed using the opacities from Fig. 2 for the models of the same four stellar types versus wavelength (horizontal axis) and geometrical height (vertical axis). Labels and overploted lines are same as in Fig. 2. The colourbar values of each star range from – max Q_{v} to max Q_{v}. Blue tones in the colourmap correspond to negative Q_{v}, i.e. cooling in the atmosphere, and red tones correspond to positive Q_{v}, i.e. heating. 
4 Opacity distribution functions
4.1 Construction of the opacity distribution function
The ODF method (Labs 1951; Strom & Kurucz 1966) was invented to allow the accurate calculation of the radiative losses due to the millions of opacity contributors in a computationally efficient way. In the 1960s the solution for the RTE was prohibitive for more than a few wavelengths in 1D model atmospheres. Despite the exponential growth in computing power since then, the same problem of efficiently and accurately accounting for the millions spectral lines is still present today when the RTE is solved in 3D timedependent simulations.
An ODF is constructed to replace detailed monochromatic opacities by their statistical means in wavelength intervals, following a procedure that preserves the accuracy of the computed bolometric radiative quantities while heavily reducing the amount of computation (see e.g. Strom & Kurucz 1966; Kurucz 1970; Vögler 2003). First, one divides the full wavelength range into steps (an example for one step, defined by λ ∈ [λ_{2}, λ_{3}], is shown at the left panel of Fig. 4). Then, all the opacities within each of the steps are sorted by increasing magnitude (see the middle panel of Fig. 4). At this point the direct mapping between the monochromatic opacity values and their wavelengths is lost. To discretise this monotonic distribution within each step, a number of substeps are defined specifying a set of dimensionless weights ω_{j}(Σ_{j} ω_{j} = 1). The length of a substep Δλ_{i,j} is calculated in terms of the length of the step Δλ_{i} and the weight ω_{j}, so that Δλ_{i,j} = ω_{j}Δλ_{i} ^{2}. Finally, the opacities are averaged computing the arithmetic mean within each substep (see the right panel of Fig. 4).
This procedure can be performed for any set of monochromatic opacities computed for a pair of thermodynamic quantities. One of them is usually the temperature T and the other is either the density ρ (as in SYNSPEC^{3}) or the gas pressure P (as in the ATLAS code, Kurucz 1970). Kurucz (1993c) created ODF tables^{4} in the described way for specific grids of (T, P). We apply the same strategy to the monochromatic opacities computed using SYNSPEC in a (T, ρ) grid, for the two monochromatic opacity tables described (see Sect. 3). The wavelength range of our SYNSPEC computation and derived ODF is shorter than the range covered by the ATLAS ODF, λ ∈ [9, 10^{7}] nm. When the ATLAS table is used to compute the bolometric RT quantities in the different model atmospheres, there is no difference between using it over its entire λ range and in the reduced range of our SYNSPEC computations. This is not surprising since the emission of cool stars peaks in spectral regions far from the extremes that we left out from the full wavelength interval of the ATLAS tables. For our ODF calculation, the wavelength range is divided into 291 steps, which have exactly the same locations and the same 12 weights for the substeps as ATLAS. The length of the steps is approximately proportional to the wavelength (see the table from Kurucz 1993c, to check the exact values); and the substeps are constructed using the weights ω_{j} = {0.1, 0.1,0.1, 0.1,0.1,0.1,0.1,0.1, 0.1,0.05,2/60,1 /60}. This kind of nonuniform weighting is taken from the ATLAS code and had been demonstrated (Cernetic et al. 2019) to perform particularly well in the visible and infrared for the solar case, being less accurate than uniform weighting in the UV (only important to compute detailed spectra). Our ODF contains the total (continuum and lines) opacity in cm^{2} g^{−1}, just as it is the case for ATLAS.
The tables from Kurucz (1993c) have been widely used in the context of MHD simulations of solar and stellar atmospheres (see, for example, Vögler 2003; Khomenko et al. 2018; Beeck et al. 2013, and many others). The atomic and molecular line lists used to compute these tables (Kurucz 1993a,b) are sufficiently complete and accurate to make these simulations closely resembling various observational diagnostics. However, apart from the molecules that are included (H_{2}, HD, MgH, NH, CH, SiH, OH, CH, C_{2}, CN, CO, SiO), many are missing (e.g. VO and TiO) which are dominant opacity sources in the spectra of M stars. SYNSPEC, however, allows for a more complete set of molecules (see Appendix A and Fig. 5) and includes an uptodate collection of atomic and molecular data (partition functions and line lists). However, to evaluate the effect of the updated data, one has to compare the two data sets carefully.
Fig. 4 Three steps in construction of the ODF. Left: schematic representation of monochromatic opacities in an interval of wavelengths divided into steps of λ_{i}. The red vertical lines mark one ODF step with λ ∈ [λ_{2}, λ_{3}]. Middle: opacities of that step sorted by magnitude. The distribution is mapped onto the unit interval w ∈ [0, 1] and the wavelength dependence is lost within the step. Right: sorted opacities discretised in ODF substeps. 
4.2 Comparison of opacity distribution function from SYNSPEC and ATLAS
While there is overall good match between the two data sets, there are also significant differences. To understand better these differences, we compare the opacity of each step and substep of the ODFs interpolated for the values of T and ρ along the four 1D model of atmospheres. Four examples representing typical differences are shown in Figs. B.1–B.4 (where black lines represent the opacity from ATLAS along the models, and the coloured lines on top of them represent the opacity from SYNSPEC).
In parallel, we also compare the bolometric Q rate computed using both ODFs. We apply the method presented in Sect. 2 to solve the RTE and compute Q_{i,j} for every step i and substep j; we then integrate in frequency following the ODF formalism:
$${Q}^{\text{ODF}}={\displaystyle \sum _{i}{Q}_{i}\equiv {\displaystyle \sum _{i}\text{\Delta}{\lambda}_{i}}{\displaystyle \sum _{j}{Q}_{i,j}{\omega}_{j}},}$$(12)
where Q_{i} is the rate integrated in wavelength in the ith step.
Comparing the values of Q computed from different data sets is not straightforward. An example of Q computed for the model of the M star using two opacity sets and their absolute difference is shown in Fig. 6. In the figure we identify the typical negative feature that corresponds to the strong cooling at the bottom of the photosphere and the relatively smaller positive one that shows the radiative heating (blanketing effect). Higher up and deeper down in the atmosphere, the values of Q can be several orders of magnitude smaller. While the radiative losses at these heights may still be significant, we focus here on the dominant components around the optical depth unity. Wherever the values of Q are small, the relative difference computed between them may be exaggerated with respect to the importance of these differences from the energy balance in the atmosphere. Moreover, as the sign of Q flips in the region of interest (the interval of heights around τ = 1) and the two Q values computed from ATLAS and SYNSPEC do not necessary flip the sign at exactly the same height, the relative difference around that height might also be misleadingly high. Therefore, instead of computing the relative differences in the entire domain, we focus on the two dominant features and measure the deviation for the cooling and heating part separately.
For each of them we compare the area A of the difference of Q computed for each of the ODF data sets normalised to the area of one of them:
$${\chi}_{\text{C}}=\frac{A{\left(\left{Q}^{\left(1\right)}{Q}^{\left(2\right)}\right\right)}_{\text{C}}}{A{\left(\left{Q}^{\left(1\right)}\right\right)}_{\text{C}}},$$(13)
$${\chi}_{\text{H}}=\frac{A{\left(\left{Q}^{\left(1\right)}{Q}^{\left(2\right)}\right\right)}_{\text{H}}}{A{\left(\left{Q}^{\left(1\right)}\right\right)}_{\text{H}}},$$(14)
where we choose as Q^{(1)} and Q^{(2)} the Q computed with the ODF from ATLAS and SYNSPEC, respectively. In Sect. 5 the same is used to measure discrepancy between the ODF and OB results. The areas are defined as $A{\left(f\left(z\right)\right)}_{C}={\displaystyle {\int}_{{z}_{\text{b}}}^{{z}_{\text{C}\to \text{H}}}f\left(z\right)\text{d}z}$ and $A{\left(f\left(z\right)\right)}_{H}={\displaystyle {\int}_{{z}_{\text{C}\to \text{H}}}^{{z}_{\text{t}}}f\left(z\right)\text{d}z}$. The height at z_{C→H} is the height where Q changes sign for the first time above the cooling component; z_{b} and z_{t} are, respectively, the highest point in the atmosphere under the surface (z < 0 in the plots) and the lowest point over the surface (z > 0) where Q < 2 × 10^{−4} min(Q) (see Fig. 6). If A (Q)_{H} < 1% A (Q)_{C}, χ_{H} is not computed (these cases are marked as χ_{H} =    % in the figures).
The radiative energy exchange terms Q_{i} computed from the four models using the same steps as in Figs. B.1–B.4 are shown in Figs. B.5–B.8. The corresponding values of the deviation measures are indicated in each panel.
When comparing the opacity from the ODF of the two data sets, we see a general match in behaviour and values in most of the visible and infrared intervals for the F, G, and K stars. An example for a step with a good match (λ ∈ [590, 600] nm) is shown in Fig. B.1. The Q_{i} values computed in the same step closely overlap for these stars with all χ_{C} and χ_{Η} between 1.4% and 3.1% (see, for example, Fig. B.5). For the nearUV, visible, and IR, the steps that include strong lines show significant differences in the opacity. Some examples are the steps that include hydrogen Hα (Fig. B.2) or NaI D1 and D2 (Fig. B.3). Similar discrepancies are found in the steps containing other strong lines like the CaII H and K lines, the G band, HI β, MgI b1, the CaII triplet, and the most intense lines from the Paschen, Brackett, and Pfund series. In the case of all the hydrogen series the opacity from ATLAS is higher than the one from SYNSPEC in the deeper part of the atmospheres, but this reverses in the upper part (see Fig. B.2). For these hydrogen lines, such opacity behaviour is present in all the substeps, although the reversal occurs deeper in the less opaque substeps. We notice a similar tendency in the G band. Another trend can be seen for the steps containing CaII lines, MgI b1 and NaI D1 and D2, for which the opacity from SYNSPEC is higher than the one from ATLAS at all the heights of the atmospheres for the most opaque substeps (see Fig. B.3). The number of affected substeps varies between the two and four most opaque, depending on which intense line falls within the step. The mismatch is affected by the number of strong lines contained in the interval, the magnitude of their opacity with respect to the continuum, and their width relative to the length of the interval. In the less opaque substeps, the two data sets give very close results.
In the steps that contain strong spectral features, the radiative energy exchange rates generally show χ_{C},χ_{H} ≲ 5–10% for the F, G, and K stars (see, for example, Figs. B.6, B.7). An exception happens for the G star for the step containing hydrogen Hα where χ_{Η} ≃ 25% (see Fig. B.6). However, in the latter example, the area of the heating feature is only 0.02 times the area of the cooling, making this discrepancy less important.
The differences between SYNSPEC and ATLAS for the strong atomic lines could be, for example, due to the wavelength resolution in the monochromatic opacity tables or the broadening parameters used in the line synthesis (e.g. Van der Waals broadening or damping constants). In the case of the hydrogen lines, the difference is explained by the treatment of the Stark broadening. For this project, we used approximate Stark broadening profiles after Hubeny et al. (1994), while in ATLAS the tables from Vidal et al. (1973) are used. Without having access to the monochromatic opacities from ATLAS it is difficult to pinpoint the exact reason of the discrepancies for the rest of the species. Although both codes are open source, their complexity makes it very difficult to find the exact source of differences by comparing the actual source code.
Apart from the spectral features discussed above, the comparison of the two opacity sets in the visible and IR for the F, G, and K stars reveals an excellent agreement. At the same time, there is a notable mismatch for the opacity in the UV between the two data sets (see an example in Fig. B.4). These discrepancies are also evident in the corresponding values of Q_{i} (see an example in Fig. B.8), for which we find differences of around χ_{C},χ_{H} ≃ 5–40% for λ ≳ 190 nm (as in the example from Fig. B.8), and around χ_{C},χ_{H} ≃ 30–80% for λ < 190 nm (in other steps in the ultraviolet).
The differences between the two data sets in the UV are likely due to the use of photoionisation crosssections from the Opacity Project^{5} and the Iron Project (Nahar 1995; Bautista 1997) in SYNSPEC (see the discussion about the missing opacity in the Sun in Allende Prieto et al. 2003a), while in the ATLAS table from 1993 the crosssections come mainly from Henry (1970) and Peach (1970).
With all the similarities found in the visible and infrared, and the nearly perfect match we see in the bolometric Q for the F, G, and K stars (Fig. 7), it may be concluded that both ODF tables are fairly similar and reproduce a close solution of the RTE.
However, both the comparison of the opacity and the computed values of Q from SYNSPEC and ATLAS reveals much larger discrepancies for the M star. The opacities from SYNSPEC are, in general, larger than the ones from ATLAS for the points at the top of the M atmosphere. The corresponding values of χ_{C} and χ_{Η} are also larger, for most of the steps in the visible and NIR they are in the range ≃10–35% (≃40–60% in the worst cases).
As previously discussed, the discrepancies between the two tables for the latertype star are expected owing to different selection of molecules and line lists. The discrepancy is prominent in the heating component of the bolometric Q (Fig. 7), for which χ_{Η} ≃ 77%, in contrast to the values χ_{C} < 3% and χ_{Η} < 10% found for the F, G, and K stars. To demonstrate the effect of certain molecules on the discrepancies between the data sets, we recompute the SYNSPEC table and Q, excluding TiO and VO from the line lists for the M star (blue curve in the bottom right panel of Fig. 7). In this case, the value of χ_{H} is considerably reduced, but it is still relatively high (≃29%). Even if the selection of molecules in the two codes is identical, differences between the opacities are expected owing to different line lists and other atomic and molecular data. For example, the use in SYNSPEC of more up to date data for the opacities, in particular for H_{2}O and TiO, from EXOMOL^{6}, or different partition functions and oscillator strengths.
Fig. 5 Oscillator strength log gf versus wavelength for all molecular spectral lines included in SYNSPEC. Notice that the transitions that mainly contribute to the opacity have typically values of log gf > −5. 
Fig. 6 Example of the radiative energy exchange rate Q versus geometrical height z for the M2V star. Top panel: the Q rate computed using the ODF from SYNSPEC data (black solid line) and the ODF from ATLAS (red solid line). The logarithm of the Rosseland optical depth is shown in the top axis. The area of the difference is shaded in grey. The zero of the geometrical height z is at the continuum optical depth at 500 nm τ_{5} = 1. The range of heights of the cooling (Q < 0) and heating (Q > 0) components are indicated in cyan and orange respectively. The heights z_{b}, z_{c→h}, and z_{t} are used as integral bounds to calculate the deviation measures χ_{C} and χ_{Η} (Eqs. (13) and (14)). The height z_{c→h} is also marked with a vertical black line. Bottom panel: absolute difference between the Q values. 
Fig. 7 Bolometric radiative energy exchange rate Q^{ODF} against the geometrical height in the four stellar models. The results using ODF from SYNSPEC (coloured lines) are compared to those using the ODF from ATLAS (thin black line). Top left of each panel: scaling factor of vertical axis. Bottom left of each panel: difference between the two ODFs computed using Eqs. (13) and (14). Top right of each panel: ratio of the heating area with respect to the cooling area for the Q for both data sets (A for ATLAS and S for SYNSPEC). In cyan in the bottom right panel for the M2V star: ODF computed using monochromatic opacity from SYNSPEC, excluding the lines from TiO and VO. The corresponding deviations for that Q curve are χ_{C} ≃ 7 and χ_{Η} ≃ 29. 
5 Opacity binning
A further approximation of the opacity computation is the OB method initially introduced by Nordlund (1982). Its implementation in the MURaM and CO^{5} BOLD codes is described by Vögler (2003) and Ludwig (1992) respectively. In this method the opacities are sorted taking into account their distribution with the discretised optical depth τ^{ref} of a representative onedimensional model atmosphere. We use the Rosseland optical depth as reference. An opacity bin l is defined by selecting a pair of lower and upper τseparators (${\tau}_{0}^{\text{ref},1},\text{\hspace{0.17em}}{\tau}_{1}^{\text{ref},\text{\hspace{0.17em}}1}$). In a variation of the method (see Magic 2014) the wavelength dependence is partially preserved by grouping opacities by both optical depth and wavelength. In that case each bin is defined by two pairs of separators, one in optical depth (τseparators) and one in wavelength (λseparators).
We start from the ODF constructed as it is described in Sect. 4. The optical depth is computed for each step i and substep j of the ODF for all geometrical heights z in a model atmosphere as
$${\tau}_{i,j}\left(z\right)={\tau}_{i,j}\left({z}_{\text{top}}\right)+{\displaystyle {\int}_{{z}_{\text{top}}}^{z}{\varkappa}_{i,j}\left(z\right)\rho \text{\hspace{0.17em}d}z},$$(15)
where ρ refers to the mass density in the model atmosphere (and not the density in the grid of the opacity table) and the height dependence in ϰ_{i,j}(z) denotes that the opacity from the ODF is interpolated into the temperatures and densities from the model atmosphere. The opacities are then grouped so that opacity of the ODF step i and substep j belong to bin l if ${\tau}_{0}^{\text{ref},1}\le {\tau}^{\text{ref}}\left({z}_{i,j}\right)<{\tau}_{1}^{\text{ref},1}$, where z_{i,j} is the height at which τ_{i,j} = 1. Additionally, when the λseparators are used, every ODF value is sorted by wavelength so that the central wavelength of an ODF step falls between a pair of λseparators. Figure 8 shows examples of such grouping of the ODF values into five bins for our four stellar models. There are two λseparators approximately dividing the wavelength axis into ultraviolet, visible, and infrared and there is one τseparator for the UV (log τ = −3) and a joint one for the visible and IR (log τ = −1). The corresponding 2D histograms of the ODF points are shown in Fig. 9. The minimal optical depth in the reference models vary with the stellar type and wavelength. We assign the top geometrical height of the atmosphere z_{i,j} = z_{top} to all ODF values with τ_{i,j}(z_{top}) > 1.
After the opacities of the ODF have been distributed into bins, for every bin l (that includes the ODF steps i = i(l) and substeps j = j(i,l) contained in the bin) and for every T, ρ of the grid of the ODF we compute the Planck function,
$${\text{B}}_{l}={\displaystyle \sum _{i\left(l\right)}\text{\Delta}{\lambda}_{i}{\text{B}}_{i}},$$(16)
its derivative,
$${\frac{\partial \text{B}}{\partial \text{T}}}_{l}={\displaystyle \sum _{i\left(l\right)}\text{\Delta}{\lambda}_{i}}\frac{\partial {\text{B}}_{i}}{\partial \text{T}},$$(17)
the Planck mean opacity,
$${\overline{\varkappa}}_{l}^{\text{Pl}}=\left({\displaystyle \sum _{i\left(l\right)}\text{\Delta}{\lambda}_{i}{\text{B}}_{i}}{\displaystyle \sum _{j\left(i,l\right)}{\omega}_{j}{\varkappa}_{i,j}}\right)/{\text{B}}_{l},$$(18)
and the Rosseland mean opacity,
$${\overline{\varkappa}}_{l}^{\text{Ro}}={\frac{\partial \text{B}}{\partial \text{T}}}_{l}/\left({\displaystyle \sum _{i\left(l\right)}\text{\Delta}{\lambda}_{i}}\frac{\partial {\text{B}}_{i}}{\partial \text{T}}{\displaystyle \sum _{j\left(i,l\right)}\frac{{\omega}_{j}}{{\varkappa}_{i,j}}}\right),$$(19)
where Δλ_{i}, ω_{j}, and ϰ_{i,j} are, respectively, the length of the wavelength step, the weight of the substep, and the opacity of the ith step and jth substep of the ODF. B_{i} is the Plank function in the middle of the wavelength step for each temperature T in the (T, ρ) grid used to compute the ODF.
Finally, following Ludwig (1992), ${\overline{\varkappa}}_{l}^{\text{PI}}$ and ${\overline{\varkappa}}_{l}^{\text{Ro}}$ are combined to get a mean value of the opacity for each bin:
$${\overline{\varkappa}}_{l}=\left({2}^{{\tau}_{l}/{\tau}_{\text{thr}}}\right){\overline{\varkappa}}_{l}^{\text{Pl}}+\left(1{2}^{{\tau}_{l}/{\tau}_{\text{thr}}}\right){\overline{\varkappa}}_{l}^{\text{Ro}},$$(20)
where the transition between thin and thick optical regimes is at the threshold optical depth τ_{thr} = 0.35 and τ_{l} is computed using pressure p in the ODF grid and a stellar surface gravity g, approximately as ${\tau}_{l}={\overline{\varkappa}}_{l}^{\text{Ro}}p/g$. With this definition, the opacity from each bin converges to the grey approximation for large enough optical depths (see the discussion about Fig. C.3).
As in the ODF case, we compute the bolometric Q rate, Q^{OB}, by solving the RTE for each bin and summing over the bins:
$${Q}^{\text{OB}}={\displaystyle \sum _{l}{Q}_{l}},$$(21)
where Q_{l} refers to partial rate of the bin l computed using Eq. (10).
We apply the described procedure to the ODF computed with SYNSPEC monochromatic opacities presented in Sect. 4. There are several important free parameters in the procedure: the model atmosphere itself, the number of bins in τ and λ, and the locations of separators between the bins. There is no obviously intuitive choice for these parameters. We designed numerical experiments to test various alternatives. As a deviation measure of the radiative energy exchange rate Q^{(2)} (computed using the binned opacity) with respect to the reference rate Q^{(1)} (computed with ODF), Vögler et al. (2004) plotted the heightdependent absolute difference of Q^{(1)} and Q^{(2)} divided by the density stratification ρ of the atmosphere. The scaling with density is important to compare the radiative rates with other energy sources and sinks in the simulations. Alternatively, Magic (2014) opted for computing max Q^{(1)} − Q^{(2)}/ max Q^{(1)}, to reduce the deviation measure to one value for a given choice of parameters. The latter choice is obviously more practical for trialanderror optimisation of the free parameters. However, in this approach the deviation is biased towards the dominant cooling component of Q. In our experiments, we take a compromise by using the deviation measures introduced in Eqs. (13) and (14). The reference values of the Q rate are those computed using ODF with SYNSPEC (Q^{(1)}) and the values that are tested (Q^{(2)}) are the ones calculated using the binned opacities. This choice of the deviation measures, allows us to automate the process of optimising the freeparameters, while it preserves separate information on both cooling and the heating component of Q. As this is a multiparameter study we break the problem down by first studying the influence of the location of the τseparators for a fixed number of bins (4) and with no bins in the wavelength direction. We then allow the number of τ bins to change. Finally, we test several cases including the λ bins.
Fig. 8 Example of the reference optical depth and wavelength of the ODF points organised into bins for the F3V, G2V, K0V, and M2V stars. As explained in the text, the Rosseland optical depth shown in the vertical axis is used as the reference τ^{ref}(z_{i,j}), where z_{i,j} is the geometrical height at which τ_{i,j} = 1 for each step and substep of the ODF. We show the same {τ, λ} binning configuration for all the stars. The colours indicate different bins. Each point corresponds to one ODF substep, all centred in the middle wavelength of their corresponding step. The last column of dots at the right of each panel corresponds to the ODF step with λ ∈ [9600, 9.5 × 10^{4}] nm, and thus, appears separated from the rest of the dots. 
Fig. 9 Histogram of the logarithm of the relative number of points from the ODF as a function of τ^{ref}(z_{i,j}) in the vertical axis and λ in the horizontal, as in Fig. 8. The density of the ODF points is the largest in the continuum (especially in the visible) and in the UV owing to the sampling of the ODF steps (see the description of the ODF tables in Sect. 4). 
Percentage of the combinations of the three separators (4 bins) with certain χ_{C} values for the four stellar types.
Percentage of the combinations of the three separators (4 bins) with certain χ_{H} values for the four stellar types.
5.1 τbin location
To understand the dependence of the Q rate on the location of the τseparators for the four stars, we first fix the number of bins to four following Nordlund (1982) and Vögler (2003). Similarly to the experiment described in Sect. 3.4 of Cernetic et al. (2019), we try all possible combinations of the τseparators from a discretised grid of optical depths for each of the stars. For the F3V star this gives us 9880 combinations selected from a grid of 40 equidistant values in log τ^{ref} ∈ [−9.0, 0.5]; for K0V and G2V, 4060 combinations (grid of 30 equidistant values in log τ^{ref} ∈ [−6.5, 0.5]); for M2V, 4060 combinations (grid of 30 equidistant values in τ^{ref} ∈ [−4.5, 0.5]). The grid is customised for each star to cover properly the relevant range in optical depths (see the vertical axis of Fig. 8). The τseparator with maximum optical depth is labelled as τ_{0}. The indices of the separators increase with height in the atmosphere (τ_{i} > τ_{i+1}).
The statistics for each star in our sample show different trends. In general, combinations of separators that give simultaneously small deviations χ^{C} and χ^{H} are easy to find for the G and K stars, while they become less frequent in the case of the F star and quite rare for the M star. The number of cases leading to the deviations below certain threshold is given in Tables 1 and 2 separately for the cooling and the heating, and in Table 3 when both components are constrained together. From Tables 1 and 2 it is clear that the cooling part is much easier to replicate accurately than the heating one. Almost any combination gives deviations less than 20% in the cooling for any of the stars, while only about half of them minimise χ_{H} below 20% for the F, G, and K and only a few for the M star. Table 3 shows how many cases we find with both χ_{C} and χ_{H} smaller than certain thresholds indicated individually for each of the stars. For the M2V star the choice of four τ bins is obviously insufficient. There are only a few combinations (less than 1%) that lead to χ_{C} < 10% and χ_{Η} < 20% and there is not a single one that reduces χ_{C} and χ_{Η} below 5 and 20% simultaneously.
To understand how the deviation measures vary with the location of the three separators, an illustrative example for the G2V star is shown in Fig. 10. The figure shows the dependence of χ_{C} and χ_{Η} on the location of the middle separator τ_{1} when the other two separators are fixed. To visualise the results, in each column of the figure the deepest separator is fixed at one of six values (log τ_{0} = 0.5, 0, −0.5, −1.2, −1.9, −2.9). Each of the curves correspond to a different fixed value of the highest separator τ_{2}. The position of the highest separator is indicated by the colour, the black curve corresponding to the smallest value of τ_{2} and the red curve to the largest one.
The figure shows that χ_{C} and χ_{H} vary significantly with the location of τ_{1} if the location of the deepest separator τ_{0} is deeper than some critical depth (lefthand columns of Fig. 10). That critical location of τ_{0} is located deeper for χ_{C} than for χ_{H} (log τ_{0} ≃ −0.5 for χ_{C} and log τ_{0} ≃ −1.2 for χ_{H}). For χ_{C} especially, if τ_{0} ≃ −0.5 the curves flatten out, meaning that the results are essentially insensitive to the location of the other two nodes.
If τ_{0} is deeper than these critical values, the deviation measures are insensitive to the location of the other two separators only if τ_{2} is relatively close to τ_{0} (the red curve in the lefthand column of Fig. 10). When τ_{0} is pushed higher up, the deviation measures are insensitive to the location of the remaining separators and their value increases (compare two righthand columns of Fig. 10).
When log τ_{0} ≥ −0.5, the shape of the curves for the χ_{C} and χ_{H} in Fig. 10 becomes more different. The cooling deviation shows one minimum around log τ_{1} ≈ −2.5, while the heating deviation has two minima: one for the locations of τ_{1} close to τ_{0} and another one that moves higher up with increasing height of τ_{0}. When log τ_{0} = −0.7, the minimum is at log τ_{1} ≈ −5, and when log τ_{0} ≥ −0.5, the minimum is at log τ_{1} ≈ −2.5. Since the location of one of the minima of χ_{H} coincides with the minimum of χ_{C}, both can be minimised using the same set of separators even when τ_{0} is as deep as log τ_{0} = 0.5.
When the location of the separator τ_{0} is set above the critical values, the sensitivity of both χ_{C} and χ_{H} to the location of τ_{1} and τ_{2} mostly disappears (see the sequence from the left to the right panel of Fig. 10). This happens for a higher up location of τ_{0} for χ_{H} than for χ_{C}^{.}
Our analysis suggests that, in the case of the G2V star, there are two strategies for placing the separators so that χ_{C} and χ_{H} remain below 4 and 15% respectively. The first strategy is to set log τ_{0} around the continuum formation height, [−0.5, 0] and log τ_{1} ∈ [−3, −2.5]. The second one is to place log τ_{0} so that it includes most of the weak and intermediate photospheric lines (log τ_{0} ∈ [−2.4, −1.9]) and log τ_{1} even higher in [−5, −3.5]. In both cases the location of log τ_{2} can be anywhere above log τ_{1}. The former case (see 2nd and 3rd column) contains the combinations that produce the minimal χ values in our experiment, although the curves are not entirely flat in such a way that the optimal solution is sensitive to the model atmosphere and to the ODF data. Alternatively, the latter case still guarantees small deviations, but it appears to be more robust as χ_{C} and χ_{H} essentially depend only on the location of the deepest bin. Since our model atmosphere does not have a chromospheric temperature rise, however, the strategy with all separators positioned relatively high may be limited to that type of models.
For the spectral types F3V and K0V the results are similar to those for G2V. For τ_{0} located deeper than a certain critical value, χ_{C} and χ_{H} vary with τ_{1} similarly to the G star, showing two minima for χ_{H}, of which one coincides with the single minimum of χ_{C}. There is a little variation of the deviation curves with τ_{2} except when it is close to the τ_{0}. For the F3V star, the optimal combinations that produce χ_{C} ≤ 8% and χ_{H} ≤ 15% are (a) log τ_{0} ∈ [−0.2, 0.5] and log τ_{1} ∈ [−2.5, −2] and (b) log τ_{0} ∈ [−2.2, −0.7] and log τ_{1} ∈ [−5.5, −4.7], with log τ_{2} again anywhere above log τ_{1} .
For the K0V star, we identified three variants that produce χ_{C} ≤ 5% and χ_{H} ≤ 15%: (a) log τ_{0} ∈ [0, 0.3] and log τ_{1} ∈ [−2, −1], (b) log τ_{0} ∈ [−0.5, −0.2] and log τ ∈ [−4, −3], and (c)log τ_{0} ∈ [−1.9, −1.7] and log τ_{1} ∈ [−5, −3], with log τ_{2} < log τ_{1}.
While the results for the M2V star presented in Fig. 11 show general similarity to the results for the other stars, there are two important differences. First, for the deep choice of log τ_{0} there is only one minimum for the χ_{H} varying with τ1 and it does not coincide anymore with the minimum for χ_{C}, but with its local maximum (see two left columns in Fig. 11). Secondly, once the τ_{0} separator is sufficiently high in the atmosphere so that the deviation measures are insensitive to changes in τ_{1}, χ_{C} flattens out as for the other stars but its value quickly increases with the height of τ_{0} (see sequence from left to right in the figure), from χ_{C} ≃ 6% when log τ_{0} = 0.2 to χ_{C} ≃ 15.5% when log τ_{0} = −0.9. Higher up, the value of χ_{C} remains larger than 12%, while χ_{H} also starts to increase once it becomes flat (after τ_{0} = −1.2). Therefore, for the M2V star it is not possible to minimise both deviation measures simultaneously with four bins. The best combination of parameters that we found gives χ_{C} ≤ 9% and χ_{H} ≤ 27% for log τ_{0} ∈ [0.2, 0.5] and log τ_{1} ∈ [−1.1, −0.6], with log τ_{2} < log τ_{1}.
Based on the analysis presented above, we selected optimal fourbins combinations as reference for other binning strategies that we described in the following sections: log τ = {−1.9, −5, −7} for the F3V star; log τ = {0, −2.5, −5} for the G2V; log τ = {−0.2, −3, −5} for the K0V; log τ = {0.3, −0.8, −2.5} for the M2V.
Fig. 10 Deviation measures χ_{C} (upper row) and χ_{Η} (lower) for the G2V star as function of the location of the separator τ_{1} when the other two separators are fixed. Columns correspond to six different fixed values of τ_{0} (values indicated at the top of each column). Each curve corresponds to a different location of the top separator τ_{2}. The colour of each curve indicates the location of τ_{2} used to compute that curve. We do not show all the possible locations for τ_{0} and τ_{2} to avoid overcrowding of the plots. The grey horizontal lines show χ_{C} = 4% and χ_{Η} = 15%. 
Percentage of the combinations of the three separators (4 bins) with certain simultaneous χ_{C} and χ_{H} values for the four stellar types.
Fig. 11 Same as in Fig. 10, but for the M2V star. The grey horizontal lines show χ_{C} = 9% and χ_{H} = 27%. 
Fig. 12 Deviation measures χ_{C} (upper panel) and χ_{H} (lower) versus the number of bins for the four stars (colours as in Fig. 1). Solid lines: deviation measures for the Q computed using 3, 4, , 11, 21, 31, 41, 51, and 101 τbins distributed equidistantly. Circles: deviations for Q computed with the Rosseland opacity. Pluses: the deviation of the optimal combinations of four τbins (see Sect. 5.1 and Table 4). 
5.2 Number of τbins
To understand the dependence of the Q rate on the number of separators, we computed Q for the four stars for different number of separators (2, 3, 5, 10, 20, 30, 40, 50, and 100)^{7} that are uniformly distributed in the relevant range for each of the stars: log τ^{ref} ∈ [−9, 0.5] for F3V, [−6.5, 0.5] for K0V, and G2V and [−4.5, 0.5] for M2V. Again, the selected optical depth ranges ensure that all the stellar models are properly represented (see the vertical axis of Fig. 8).
In Fig. 12 χ_{C} and χ_{H} are shown as a function of the number of bins. For F, G, and K  type stars, there is very little variation in χ_{C} The trend is slightly different for each of the stars individually. While χ_{C} slowly decreases for the K star with an increasing number of bins, it increases for G2V, and it first increases and then decreases for the F3V – type star. In all three cases the value of χ_{C} saturates at around 20 bins. The M star shows a significant improvement in χ_{C} when the number of bins increases from 4 to 21. For more bins it saturates at ≈5%. Although there is more variation in χ_{Η} than in χ_{C}, the overall trend is similar and, again, the results saturate for more than ≈20 bins. Vögler et al. (2004) tested the difference between the Q rate computed with the binned opacity and with the ODF for an increasing number of bins (up to 60). Similarly to our results, they found that the solutions for Q converge towards a limiting value that is different from Q_{ODF}.
Wray et al. (2006) investigated the performance of the OB and ODF methods in the context of shockgenerated radiation during simulated reentry of Apollo AS501 vehicle into the Earth’s atmosphere. Although the physics of the experiment and implementation of the OB are different (e.g. Planckian mean for ${\overline{\varkappa}}_{l}$, ϰsorting instead of τsorting), they also found that the results of the radiation computation quickly converge with increasing number of bins (reaching saturation at around 10 bins). However, note that their conclusions about the ODF solution are not directly translatable to the radiative transfer in the stellar atmospheres owing to different definitions of ODF.
Therefore, the saturation values of χ_{C} and χ_{H} are the minimal deviations that cannot be eliminated by increasing the number of the τseparators. Depending on the stellar type, the minimal χ_{C} values span between 2 and 7%, and the minimal χ_{H} span between ≈10 and ≈30%. These minimal deviations are the smallest for the K and the largest for the F star. However, it is possible to find certain combinations of the separators that minimise the deviations particularly well. For example, the χ_{H} value for the F3V star is reduced to only a few per cent when three separators are set to log τ = {−1.8, −4.3, −6.6}. This is consistent with locations that produce simultaneous minimum χ_{C} and χ_{H} deviations, as discussed in Sect. 5.1. The deviations with the optimal combinations found using four bins for each of the stars in Sect. 5.1 are shown as coloured plus markers in Fig. 12. These deviations are either similar to or smaller than the saturation values when the bins are distributed uniformly. We have also computed the Q rates using the Rosseland opacity mean applied to the whole spectrum (Eq. (C.1)). The χ differences for this grey solution are shown as coloured circles in Fig. 12. Their values are larger or similar to the saturation ones, but always larger than the solutions with four bins from Sect. 5.1.
We experimented with several other strategies for distributing separators automatically. They all confirm the conclusion that careful distribution of the separators is more important than their number.
Fig. 13 Example of the binning used to compute the Q values for the G2V star, using different numbers of {τ, λ}bins. Similarly to Fig. 8, each grey dot corresponds to one ODF substep. For this star, the τseparators in the UV in all the experiments are fixed at log τ = {−1.8, −4.2}. Columns, from left to right, show the cases for N = 6–9 bins. The red numbers show indices of the bins. The case for N = 6 has two λseparators at 400 and 1600 nm. For more than 6 bins, bins 1 to 5 are kept identical, but more λ separators are added between 400 and 9750 nm. A similar binning setup is used for the other stars with modified τseparators in the UV and the λ separator between the short and long wavelengths for N > 6. 
5.3 {τ, λ} binning
The OB method using only τseparators does not reproduce the exact solution in the limit of infinite number of bins because the contributions from different parts of the spectrum, each with its own height variation, are mixed in any given bin (Vögler et al. 2004, see discussion in Appendix C). One way to improve the approximation is to preserve some of the wavelength dependence in the binning procedure by grouping the opacity data points with respect to their spectral regions prior to the τbinning. Therefore, λseparators between these regions are introduced in addition to the rseparators. Such bidimensional distribution of the opacities is nonuniform in (τ, λ) parameter space: in every wavelength interval the τseparators may be specified independently (both their number and location) to account for the height distribution of the opacities in that interval. Ludwig (1992) introduced the {τ, λ}binning, showing that splitting the continuum at the Balmer jump improved significantly the results for Vega, due to the systematically different behaviour of Q for wavelengths shortwards and longwards of the jump. Vögler (2003) explored the possibility of splitting the least opaque bin into two subgroups for long and short wavelengths and concluded that the improvement is probably not worth the extra computational effort. Collet et al. (2018) also attempted to improve the OB procedure by partly accounting for the lost wavelength dependence. Magic (2014) treated the UV and the visible + IR opacities separately. Each group is sorted with respect to its own set of τseparators. In addition, the least opaque bin of the visible + IR opacities is further split with a set λseparators.
We computed the Q values using a {τ, λ}binning approach similar to Magic (2014, see his Fig. 2.6) for different total number of bins, N = 6, 7, 8, 9, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, and 100. The distribution of the bins for the G2V star for N = 6, 7, 8 and 9 is shown in Fig. 13. For any N, the spectrum is first divided into the short and long wavelength parts by placing one λseparator at λ_{0}. The value of λ_{0} is 400 nm for the F and G stars (effectively separating UV from visible), 500 nm for the K star, and 600 nm for the M star. The shortwavelength opacities are then grouped into three τbins, and the longwavelength ones into two. The least opaque τbin of the longwavelength opacities is then split by one or more additional λseparators. For N = 6 the λseparator is placed at 1600 nm. For N > 6, the λseparators of the least opaque τbin are placed equidistantly in log λ in the range [λ_{0}, 9570] nm. The τseparators in the shortwavelength region are customised for each stellar type: −2.6 and −5.8 for F3V, −1.8 and −4.2 for G2V, −1.8 and −4.0 for K0V, and −1.2 and −2.8 for M2V. The locations of the separators are chosen optimally to sample different features present for each of the stars when the ODF points are plotted in the τλ plane (cf. Figs. 8 and 9).
We computed the Q rate for every combination of {τ, λ}bins. The corresponding deviations χ_{C} and χ_{H} are shown in Fig. 14. The behaviour of χ_{C} with the number of bins is generally similar to that found in Fig. 12. For four stars the χ_{C} deviation decreases with increasing number of bins until it saturates after around 10–20 bins. The saturation values of χ_{C} are lower than the ones found in Sect. 5.2 (cf. Table 4). For less than 15 bins (Δ log λ > 0.1) χ_{C} remains nearly constant for G2V and K0V, but increases significantly for F3V and M2V. The results for M2V show two peaks for N = 7 and 10, and a minimum between them for N = 8.
This strong variation in χ_{C} is a good illustration of how sensitive the Q rate is on how the ODF data points are distributed between the bins. Figure 13 shows the distribution of the bins with N = 7 and N = 9 (2nd and 4th panel). Bins 1 to 5 are identical in both cases. Bins 6 and 7 for N = 7 are split into two bins each for N = 9. Bin 7 (for N = 7) includes most of the visible spectrum, it coincides with the central part of the Planck function and thus it dominates the radiative losses. The rate Q computed from that bin is compared to the total Q in the lefthand panel of Fig. 15. In the righthand panel it is compared with Q values computed from the bins 8 and 9 (for N = 9) and with their sum. Individual Q profiles for the bins 8 and 9 (for N = 9) differ significantly in shape, amplitude and location of the maximal cooling. Because of the nonlinearity of RTE, their sum is not equal to Q computed when the opacity of the two bins is merged in the bin 7 for N = 7. In terms of our deviation measure χ_{C}, this leads to a drop from ≈14% for N = 7 to ≈3% for N = 9 (top panel of Fig. 14). The finer the distribution of the visible opacities in wavelength is, the more stable are the results leading to the saturation of χ_{C}.
For the M2V star the χ_{H} saturation value is about the same as in Fig. 12, but for the F, G, and K stars χ_{H} saturates at significantly larger values owing to an excess of heating relative to the ODF solution.
In another experiment, we tested a variant of the previous method were bins are distributed in a checkerboard pattern in the (τ, λ) plane. The number of τseparators is varied between 1 and 16 and the number of λseparators is fixed at 5 (λ = 381, 726, 1383, 2636, 5023 nm). The τseparators for all the stars are distributed equidistantly at log τ^{ref} ∈ [−3, 0.5]. The minimum number of bins in this experiment is 12, the maximum 102. The results are shown in Fig. 16. For a small number of bins χ_{H} is larger than the deviation measure we obtain with the optimised solutions from Sect. 5.1. With 12 bins the values χ_{H} lie between 20 and 30%, similar to those in Fig. 12 with 11 bins distributed only in τ. With finer sampling in τ, the value of χ_{H} is quickly reduced below 10% and saturated at that level for all four stars.
From this analysis it follows that introduction of the λ separators may reduce the minimum deviation of the OB in particular cases when it is combined with carefully selected separators in τ. However, this improvement is relatively small in comparison with the basic method where sorting is done only by optical depths (cf. Figs. 16 and 12). Table 4 summarises the deviation measures χ_{C} and χ_{H} computed in the experiments described in this section. The values shown for the experiments described in Sects. 5.2 and 5.3 correspond to the minimum number of bins at which the deviation measures approach their saturation values. It is clear that the 4bin approach with carefully tuned locations of the τseparators (Sect. 5.1) performs similarly or better than the other approaches for all four stars. The only exception is the M star, where taking into account the wavelength distribution simultaneously with distribution in τ for all λ intervals is an obvious advantage (Fig. 14). The strategy proposed by Magic et al. (2013) reproduces the cooling component better than other strategies, but it produces large errors in the heating component.
Finally, in Fig. 17 Q and Q/ρ are compared for the four stars and five opacity approaches (ODF, grey opacity, and three OB setups). The deeper part of the atmosphere where mainly cooling occurs is shown in the two left columns; the higher layers where the heating occurs are shown in the two right columns. On one hand, the cooling can be fitted either using Q or Q/ρ, since the shape of Q is practically not modified by the relatively slowly changing density in subsurface layers. On the other hand, the density scale height diminishes faster in higher layers, thus, dividing by the rapidly decreasing density enhances the differences we see in Q. For the part of the atmosphere studied in the present paper (top of the convection zone and lower photosphere), Q and Q/ρ show similar information regarding the quality of the computed Q^{OB} in respect to Q^{ODF}. Some differences are expected when choosing Q/ρ, but the general behaviours described in this work should be the same. If one is interested in layers closer to the chromosphere and beyond, Q/ρ would be the quantity to optimise.
Fig. 14 Deviation measures χ_{C} (upper panel) and χ_{Η} (lower) versus the number of bins, for the four stars (colours as in Fig. 1), for the case of {τ, λ)binning approach similar to Magic (2014). Solid lines: deviation measures for the Q computed using 6, 7, 8, 9, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, and 100 bins distributed as in Fig. 13. Circles: deviations for Q computed with the Rosseland opacity. Crosses: deviations for Q computed with only two τseparators located at the heights of the bins in the UV (no λseparators). Pluses: the deviation of the optimal combinations of four τbins (see Sect. 5.1 and Table 4). 
Deviation measures χ_{C} and χ_{H} for different binning strategies described in Sects. 5.1–5.3, for the four stars F3V, G2V, K0V, and M2V (columns are named after the corresponding sections).
Fig. 15 Comparison of the radiative energy exchange rate Q from bin 7 for N = 7 and bins 8 and 9 for N = 9 corresponding to the cases shown in Fig. 14 for the M2V star. Left panel: Q_{7} dominates the total Q^{OB} = Σ_{l} Q_{l} for N = 7. Right panel: comparison of Q_{8} and Q_{9} for N = 9 against Q_{7} for N = 7. 
Fig. 16 Deviation measures χ_{C} (upper panel) and χ_{Η} (lower) against the number of bins for the four stars, computed for the strategy with five fixed λseparators at λ ∈ (381, 726, 1383, 2636, 5023) nm and with varying number of τseparators (1, 2, 3, 5, 6, 8, 10, 12, 14, and 16) and therefore varying total number of bins (12, 18, 24, 36, 42, 54, 66, 78, 90, and 102 bins). Solid lines: deviation measures for Q computed using a different number of {τ, λ}bins. Circles: deviations for Q computed with the Rosseland opacity. Pluses: deviations of the optimal combinations of four τbins (see Sect. 5.1 and Table 4). 
6 Conclusions
This work is based on detailed monochromatic opacity tables computed with the SYNSPEC (Hubeny & Lanz 2011) code using complete uptodate sets of atomic and molecular line lists. The ODF is calculated using these opacities and compared with the ODF from the ATLAS code (Kurucz 1993c). Although there are some systematic differences between the two data sets, a good match is found for the F3V, G2V, and K0V stellar spectral types when the radiative energy exchange rate Q is computed using each of the data sets. For the M2V – type star the comparison of SYNSPEC versus ATLAS ODFs reveals the importance of the molecules in the spectra of these stars. The lack of some significant opacity contributors such as VO and TiO in the ATLAS line lists produces less heating in the radiative energy exchange rate than in the SYNSPEC results.
In Sect. 5 the OB method (Nordlund 1982) is applied to the SYNSPEC ODF. This method enabled rapid progress in realistic 3D MHD simulations of the nearsurface convection as it reduced the number of required radiative transfer computations by several orders of magnitude while preserving the essence of the nongrey nature of the radiative transfer solution. However, the method involves several free parameters that cannot be intuitively guessed. We systematically tested several strategies for the opacity binning: positioning of the four bins (Sect. 5.1), the number of τ bins (Sect. 5.2) and simultaneously binning in optical depth and wavelength (Sect. 5.3). None of the strategies tested converge to the ODF solution in the limit of a large number of bins since the binning mixes parts of the spectrum with different height profiles of the opacity (Appendix C). It is possible to find combinations of small number of bins that significantly reduce the deviations. However, these combinations are modeldependent and with relatively narrow intervals in which the bins can be located. Strategies that include preservation of the wavelength dependence through combined τ and λbinning may improve the deviation measures, but with a significant increase in computing overheads. For example, even in the most extreme case of the M star, to reduce the already small χ values by a factor of three (see the bottom row of Table 4), requires an increase in the number of bins, and thus the computing cost by a factor of four or five. In most applications such accuracy is not needed as other approximations made in the RTE solution and the MHD modeling contribute more to the total error budget.
The results presented in this study are limited to the solar chemical composition. In the stars with different metallicities both the atmospheric structure and the relative importance of the various opacity contributors will change. Therefore, our results should not be blindly applied to these stars. Nevertheless, the trends identified in our analysis may be used as guidelines for tuning the details of the OB procedure for the stars with nonsolar metallicities.
An open question left in this work is the importance of the measured deviations χ_{C} and χ_{H} in the final structure of the atmosphere in the context of simulations. To answer this question, one may start by extending the study for the Sun in Sect. 3.4 of Vögler et al. (2004) for the case of other stellar spectral types in time dependent simulations.
It should also be emphasised that the optimal binning approximation selected based on the 1D models will not necessarily remain optimal when implemented in 3D simulations. However, the tuning of the binning strategy is computationally too expensive to be routinely done in 3D. In the followup of this study we shall explore how different binning strategies perform in realistic 3D atomospheres simulated with the MANCHA code (Khomenko et al. 2017, 2018), especially in the case of the M – type star.
Fig. 17 Comparison of Q and Q/ρ for the four stars (each row) and five opacity approaches: ODF (black solid line), grey opacity computed with Eq. (C.1; grey solid line), and the binned opacities for four bins (red solid line), 15 bins (yellow solid line), and 24 bins (black dotted line) corresponding to Cols. 5.1, 5.3a and 5.3b in Table 4. Left column: Q for the deeper part of the atmosphere where mainly cooling occurs. Second column: Q/ρ for the same part of the atmosphere. Third column: Q for the higher layers where the heating occurs. Right column: Q/ρ for the same part of the atmosphere. 
Acknowledgements
This work was supported by the European Research Council through the Consolidator Grant ERC–2017–CoG–771310–PI2FA and by Spanish Ministry of Science through the grant PID2021–127487NB–I00. A.P.G. acknowledges support from the Agencia Estatal de Investigación (AEI) of the Ministerio de Ciencia, Innovación y Universidades (MCIU) and the European Social Fund (ESF) under grant with reference PRE2018–086567. C.A.P. is thankful for funding from the Spanish goverment through grants AYA2014–56359–P, AYA2017–86389–P and PID2020–117493GB–100. We thank the referee for the detailed and careful reading of the manuscript and for asking the questions that lead to Appendix C and Fig. 17. We thank Hans Günter Ludwig for reading the manuscript, the interesting discussions about the OB method and for suggesting the analysis in Fig. C.4. We are thankful to Alexander Shapiro for reading the manuscript and giving his feedback. We are thankful to Terry Mahoney for the language editing. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.
Appendix A Content of the SYNSPEC monochromatic opacity
In SYNSPEC each of the transitions taken into account to compute the opacity can be specified by the user. Only explicit ions (in our computations, H, H^{+}, H^{−}, He, He^{+}, He^{++} C, C^{+}, C^{++}, N, N^{+}, N^{++}, O, O^{+}, O^{++}, Na, Na^{+}, Na^{++}, Mg, Mg^{+}, Mg^{++}, Al, Al^{+}, Al^{++}, Si, Si^{+}, Si^{++}, Ca, Ca^{+}, Ca^{++}, Fe, Fe^{+} and Fe^{++}) are considered in the continuum opacities solved in LTE (see sections 2.4 and 2.5 of Hubeny & Lanz 2017a as well as Hubeny et al. 2021 for more details). We use the data included in the synple repository (version 1.2), as described in the following paragraphs. The contributions involving molecules and any of the CollisionInduced Absorptions are only taken into account below certain temperature, in this work chosen to be 8000K (i.e. when molecules are considered in the EOS).
The freefree crosssections for H and He^{+} are computed hydrogenically with approximated freefree Gaunt factor (see chapter 7.1 of Hubeny & Mihalas 2014). For He, C, N, O, Na, Mg, Al, Si, Ca and Fe the crosssections are computed again hydrogenically, but with Gaunt factor of 1. Other freefree contributions are also taken into account:
H^{−} from Bell & Berrington (1987).
${\text{H}}_{2}^{+}$ from Gingerich (1969).
${\text{H}}_{2}^{}$ from Bell (1980).
He^{−} from the polynomial fit of Carbon et al. (1969) to the data of John (1968).
The boundfree transitions are mainly computed using hydrogenic crosssection, data from the Opacity Project TOPbase^{8} and the Iron Project ^{9}:
H^{−}: the crosssection is computed hydrogenically with Gaunt factor of 1 (see chapter 7.1 of Hubeny & Mihalas 2014).
${\text{H}}_{2}^{+}$: the crosssection is computed using the expression from Gingerich (1969).
H: the crosssection for the first nine transitions are computed hydrogenically with Gaunt factor of 1.
He: the crosssections for the first 14 transitions are computed using cubic fits from Fernley et al. (1987) from the Opacity Project data, taking into account the multiplicity of every transition (if these are triples or singlets).
He^{+} : the crosssection for 14 transitions is computed hydrogenically with exact Gaunt factors.
C: the crosssections for 104 transitions are interpolated from TOPbase data.
C^{+} : the crosssections for 40 transitions are interpolated from TOPbase data.
N: the crosssections for 89 transitions are interpolated from TOPbase data.
N^{+}: the crosssections for the 51 transitions are interpolated from TOPbase data.
O: the crosssections for 54 transitions are interpolated from TOPbase data.
O^{+}: the crosssections for 74 transitions are interpolated from TOPbase data.
Na: the crosssections for 32 transitions are interpolated from TOPbase data.
Na^{+}: the crosssections for eight transitions are interpolated from TOPbase data.
Mg: the crosssections for 71 transitions are interpolated from TOPbase data.
Mg^{+}: the crosssections for the 31 transitions are interpolated from TOPbase data.
Al: the crosssections for 33 transitions are interpolated from TOPbase data.
Al^{+} : the crosssections for 81 transitions are interpolated from TOPbase data.
Si: the crosssections for 57 transitions are interpolated from TOPbase data.
Si^{+}: the crosssections for 46 transitions are interpolated from TOPbase data.
Ca: the crosssections for 79 transitions are interpolated from TOPbase data.
Ca^{+}: the crosssections for 32 transitions are interpolated from TOPbase data.
Fe: the crosssections for 49 transitions are obtained from the Iron Project data.
Fe^{+}: the crosssections for 41 transitions are obtained from the Iron Project data.
Other contributions to the continuum and bands are also included:
H, He and H_{2} Rayleigh scattering from Dalgarno (1962).
Thomson scattering: ϰ = 6.65 × 10^{−25}n_{e}/ρ (cm^{2}g^{−1})
CollisionInduced Absorption of H_{2}H_{2} (Borysow et al. 2001), H_{2}He (Jørgensen et al. 2000), H_{2} H (Gustafsson & Frommhold 2003) and HHe (Gustafsson & Frommhold 2001).
CH, OH continuous absorption from Kurucz et al. (1987).
Boundbound transitions are computed using the data from different line lists, from which the code selects which lines may actually contribute (e.g. with a threshold in the linetocontinuum opacity ratio). These are mainly from Kurucz line lists ^{10}, updated by data from the National Institute of Standards and Technology (NIST) when available. SYNSPEC also uses data from EXOMOL^{11} for some of the molecular transitions. In the case of the molecules, around a total of 2 × 10^{7} transitions for H_{2}, CH, NH, OH, NaH, MgH, SiH, CaH, CrH, FeH, C_{2}, CN, CO, MgO, AlO, SiO, CaO, VO, H_{2}O and TiO are included (see Fig. 5). In the case of the atoms, SYNSPEC includes a line list with around 2 × 10^{6} transitions for:
Neutral H, As, Se, Rb, Sb, Te, Cs, Pt, Au, Tl, and Bi.
Neutral and first ion of He, Li, Ga, Ge, Sr, Y, Tc, Pd, Ag, Cd, In, Sn, Ba, La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Tm, Yb, Lu, Hf, Ta, W, Re, Os, Ir, Hg, Pb, Th, and U.
Neutral and first two ions of Be, Nb, Rh, Zr, and Ru.
Neutral and first three ions of B, C, and Mo.
Neutral and first five ions of N, O, F, Ne, Na, Mg, Al, Si, P, S, Cl, Ar, and K.
Neutral and first seven ions of Zn.
Neutral and first eight ions of Ca, Sc, Ti, V, Cr, Mn, and Co.
Neutral and first nine ions of Fe, Ni, and Cu.
Appendix B Details of the comparison of the opacity distribution function from SYNSPEC and ATLAS
Figures in this appendix show examples of typical differences between the ODFs from SYNSPEC and ATLAS in the opacity and in the Q rates for the four stellar models (see Sect. 4.2).
Fig. B.1 ODF opacity against the temperature of the four stellar atmospheres (different colour lines, corresponding to a different spectral type as in Fig. 1), and for 12 substeps (different panels) of the ODF step defined by the λ ∈ [590, 600] nm range: SYNPEC (thick coloured lines) and ATLAS (thin black lines). The number at the left bottom corner of each panel indicates the substep in the order of increasing opacities. This ODF step is an example of a generally good match between SYNPEC and ATLAS that is found in the majority of ODF steps in the visible and IR. The only systematic difference is present in all substeps at the top part of the M model. 
Fig. B.2 Same as Fig. B.1, but for the step defined by λ ∈ [650, 660] nm containing hydrogen Ha. This example shows that there are significant differences between the two data sets in the strong hydrogen line. The differences are present in all substeps for the M star, but also in the more opaque substeps (rows three and four) for the three hotter stellar types. 
Fig. B.3 Same as Fig. B.1, but for the step defined by λ ∈ [580, 590] nm containing NaI D1 and D2 line. 
Fig. B.4 Same as Fig. B.1, but for the step defined by λ ∈ [193, 197.5] nm. The differences in the UV between the two tables are large in all the substeps and for all stellar models. 
Fig. B.5 Radiative energy exchange rate Q_{i} computed from the four model atmospheres using the opacity shown in Fig. B.1 and integrated over the interval λ ∈ [590, 600] nm. The values on the vertical axis are scaled with the factor indicated in the top left of each panel. Bottom left of each panel: deviation between the two ODFs computed using Eqs. 13, 14. Top right of each panel: ratio of the heating area with respect to the cooling area for the Q_{i} for both data sets (A for ATLAS and S for SYNSPEC). 
Fig. B.6 Same as Fig. B.5, but for the λ ∈ [650, 660] nm interval that contains the hydrogen Hα line. Computed using the opacity shown in Fig. B.2. 
Fig. B.7 Same as Fig. B.5, but for λ ∈ [580, 590] nm, interval that contains the NaI D1 and D2 line. Computed using the opacity shown in Fig. B.3. 
Fig. B.8 Same as Fig. B.5, but for λ ∈ [193, 197.5] nm, computed using the opacity shown in Fig. B.4. 
Appendix C Saturation of the opacity binning method and comparison against grey opacity
The OB method is not expected to converge to the ODF solution for a large number of τbins (Vögler et al. 2004). To explain why this is the case, we use the checkerboard distribution of the bins introduced in Fig. C.1, similar to that used in Sect. 5.3, but now with four bins in τ and six in λ. The number of nonempty τbins changes depending on the λbin. For example, in the case of the F3V star there are four τbins for λ < 300 nm, while there are two τbins for λ ∈ [1000, 2000] nm (Fig. C.1). For each star, in Fig. C.2 we compare the binned opacity of pairs of wavelength ranges selected to better illustrate this behaviour. We compare λranges that have the same number and location of nonempty τbins. For each star and any of the τbins, the opacity from one of the two λranges vary with height in a different way than the opacity from the other. Even for arbitrarily increased number of bins, binning only in τ mixes wavelength regions which opacities have different height profiles (this is illustrated in figure 4.13 from Vögler 2003, too). Although this effect may be reduced, it is present even when splitting in λ, since each contributor to the opacity may have its own variation with height.
Figure C.3 shows how the opacity computed with the OB method converges to the Rosseland opacity in each τbin, ${\overline{\varkappa}}_{l}={\overline{\varkappa}}_{l}^{\text{Ro}}$ (dotted black lines), in the limit of large enough optical depths. We show only the case for the F3V star and the range for λ < 300 nm to avoid overcrowding the plot. However, the conclusions are valid for any λrange and star. Another example of this convergence can be found in figure 4.12 from Vögler (2003) for the case of the Sun. Figure C.3 also shows the Rosseland opacity mean for the whole spectrum in black dashed line. The harmonic mean used to compute this averaged opacity,
$${\overline{\varkappa}}_{\text{Ro}}={\displaystyle \int \frac{\partial \text{B}}{\partial \text{T}}d\lambda /\left({\displaystyle \int \frac{\partial \text{B}}{\partial \text{T}}\frac{1}{{\varkappa}_{\lambda}}d\lambda}\right),}$$(C.1)
gives more weight to the continuum than to the line cores, making that the Rosseland mean is close to the opacity of the least opaque bins.
Finally, Figure C.4 shows that the combination of the opacity bins in the Rosseland fashion,
$$\sum _{l}{\frac{\partial B}{\partial T}}_{l}/\left({\displaystyle \sum _{l}{\frac{\partial B}{\partial T}}_{l}\frac{1}{{\overline{\varkappa}}_{l}}}\right)},$$(C.2)
yields the correct value for the Rosseland mean for deep enough optical depths. For these optical depths (see Fig. C.3 and Eq. 20)
$${\overline{\varkappa}}_{l}\approx {\overline{\varkappa}}_{l}^{\text{Ro}}={\frac{\partial \text{B}}{\partial \text{T}}}_{l}/\left({\displaystyle \sum _{i\left(l\right)}\text{\Delta}{\lambda}_{i}}\frac{\partial {\text{B}}_{i}}{\partial \text{T}}{\displaystyle \sum _{j\left(i,l\right)}\frac{{\omega}_{j}}{{\varkappa}_{i,j}}}\right),$$(C.3)
thus, the combination of the bins following Eq. C.2 gives
$$\begin{array}{l}{\displaystyle \sum _{l}{\frac{\partial B}{\partial T}}_{l}/\left({\displaystyle \sum _{l}{\frac{\partial B}{\partial T}}_{l}\frac{1}{{\overline{\varkappa}}_{l}}}\right)}=\hfill \\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\displaystyle \sum _{l}{\displaystyle \sum _{i\left(l\right)}\text{\Delta}{\lambda}_{i}}}\frac{\partial {\text{B}}_{i}}{\partial \text{T}}/\left({\displaystyle \sum _{l}{\displaystyle \sum _{i\left(l\right)}\frac{\partial {\text{B}}_{i}}{\partial \text{T}}}{\displaystyle \sum _{j\left(i,l\right)}\frac{{\omega}_{j}}{{\varkappa}_{i,j}}}}\right).\hfill \end{array}$$(C.4)
This equation is similar to Eq. C.1, but with the formalism of the ODF to compute the integrals and with the additional sums over the bins. It can be written in a more general way as
$$\begin{array}{l}{\displaystyle \sum _{l}{\frac{\partial B}{\partial T}}_{l}/\left({\displaystyle \sum _{l}{\frac{\partial B}{\partial T}}_{l}\frac{1}{{\overline{\varkappa}}_{l}}}\right)}\equiv \hfill \\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\displaystyle \sum _{l}{\displaystyle {\int}_{\lambda \left(l\right)}\frac{\partial {\text{B}}_{\lambda}}{\partial \text{T}}d\lambda}}/\left({\displaystyle \sum _{l}{\displaystyle {\int}_{\lambda \left(l\right)}\frac{\partial {\text{B}}_{\lambda}}{\partial \text{T}}\frac{1}{{\varkappa}_{\lambda \left(l\right)}}}}d\lambda \right),\hfill \end{array}$$(C.5)
where λ(l) corresponds to the wavelengths that enter the lth bin. Binning in τ, λ, or both always is equivalent to group only separating in λ, since each opacity point corresponds to a different wavelength^{12}. This property as well as understanding the integrals of the opacity in terms of Riemann or Lebesque integration explains the result shown in Fig. C.4. Similarly, the Planck mean should be recovered for higher layers.
Fig. C.1 Distribution of the opacity bins in optical depth and wavelength for the four stars. Similarly to Fig. 8, each coloured dot corresponds to one ODF substep. The wavelength separators are identical for all stars (λ < 300, ∈ [300, 500], ∈ [500, 700], ∈ [700, 1000], ∈ [1000 – 2000], and > 2000 nm). All values within one wavelength range are shown in the same colour. The τseparators are selected from the optimal combinations found in Sect. 5.1 for each star: log τ = {−1.9, −5, −7} for the F3V star; log τ = {0, −2.5, −5} for the G2V; log τ = {−0.2, −3, −5} for the K0V; log τ = {0.3, −0.8, −2.5} for the M2V. 
Fig. C.2 Opacity ${\overline{\varkappa}}_{l}$ of each τbin in pairs of λranges for the four stars (spectral type at the top left of each panel) versus the geometrical height (bottom axis of each panel) and optical depth (top axis). The λranges are selected to illustrate the different variations of the opacity with height depending on the spectral region. The colours of the curves correspond to the λranges in Fig. C.1, with the values of each range in the legend. 
Fig. C.3 Opacity ${\overline{\varkappa}}_{l}$ for each τbin in λ < 300 nm computed with Eq. 20 (solid blue curves) compared to the opacity ${\overline{\varkappa}}_{l}={\overline{\varkappa}}_{l}^{\text{Ro}}$ (dotted black curves) for the F3V star. The rest of the coloured curves correspond to the opacity for the least opaque bin for each wavelength range. Dashed black line: Rosseland opacity mean computed with Eq. C.1. The bottom and top axis show the geometrical height and optical depth, respectively. 
Fig. C.4 Rosseland opacity mean computed with Eq. C.1 (black dashed line) and Eq. C.2 (red thin solid line) for the F3V star. The bottom and top axis show the geometrical height and optical depth, respectively. 
References
 Allende Prieto, C., Hubeny, I., & Lambert, D. L. 2003a, ApJ, 591, 1192 [NASA ADS] [CrossRef] [Google Scholar]
 Allende Prieto, C., Lambert, D. L., Hubeny, I., & Lanz, T. 2003b, ApJS, 147, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Allende Prieto, C., Hubeny, I., Lanz, T., & Osorio, Y. 2022, Synple, https://github.com/callendeprieto/synple/blob/master/docs/synple.pdf [Google Scholar]
 Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
 Bautista, M. A. 1997, A&AS, 122, 167 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beeck, B., Cameron, R. H., Reiners, A., & Schüssler, M. 2013, A&A, 558, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bell, K. L. 1980, J. Phys. B: Atom. Mol. Phys., 13, 1859 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, K. L., & Berrington, K. A. 1987, J. Phys. B: Atom. Mol. Phys., 20, 801 [NASA ADS] [CrossRef] [Google Scholar]
 Borysow, A., Jorgensen, U. G., & Fu, Y. 2001, J. Quant. Spec. Radiat. Transf., 68, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233 [NASA ADS] [Google Scholar]
 Caffau, E., Ludwig, H. G., Steffen, M., et al. 2008, A&A, 488, 1031 [CrossRef] [EDP Sciences] [Google Scholar]
 Carbon, D., Gingerich, O. J., & Latham, D. W. 1969, in LowLuminosity Stars, ed. S. S. Kumar (USA: Gordon and Breach), 435 [Google Scholar]
 Castelli, F., & Kurucz, R. L. 2003, IAU Symp., 210, A20 [Google Scholar]
 Cernetic, M., Shapiro, A. I., Witzke, V., et al. 2019, A&A, 627, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Collet, R., Nordlund, Å., Asplund, M., Hayek, W., & Trampedach, R. 2018, MNRAS, 475, 3369 [NASA ADS] [CrossRef] [Google Scholar]
 Dalgarno, A. 1962, Geophys. Corp. of America, Tech. Rep., 62 [Google Scholar]
 Felipe, T., Khomenko, E., & Collados, M. 2010, ApJ, 719, 357 [Google Scholar]
 Fernley, J. A., Seaton, M. J., & Taylor, K. T. 1987, J. Phys. B: Atom. Mol. Phys., 20, 6457 [NASA ADS] [CrossRef] [Google Scholar]
 Gingerich, O. 1969, Theory and Observation of Normal Stellar Atmospheres (Cambridge: M.I.T. Press), 377 [Google Scholar]
 Gingerich, O., Noyes, R. W., Kalkofen, W., & Cuny, Y. 1971, Sol. Phys., 18, 347 [Google Scholar]
 Gustafsson, M., & Frommhold, L. 2001, ApJ, 546, 1168 [NASA ADS] [CrossRef] [Google Scholar]
 Gustafsson, M., & Frommhold, L. 2003, A&A, 400, 1161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 42, 407 [NASA ADS] [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayek, W. S. 2010, PhD thesis, Australian National University, Canberra, Australia [Google Scholar]
 Henry, R. J. W. 1970, ApJ, 161, 1153 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., & Lanz, T. 2011, Astrophysics Source Code Library [record ascl:1109.022] [Google Scholar]
 Hubeny, I., & Lanz, T. 2017a, arXiv eprints [arXiv:1706.01935] [Google Scholar]
 Hubeny, I., & Lanz, T. 2017b, arXiv eprints [arXiv:1706.01937] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton: Princeton University Press) [Google Scholar]
 Hubeny, I., Hummer, D. G., & Lanz, T. 1994, A&A, 282, 151 [NASA ADS] [Google Scholar]
 Hubeny, I., Allende Prieto, C., Osorio, Y., & Lanz, T. 2021, arXiv eprints [arXiv:2104.02829] [Google Scholar]
 John, T. L. 1968, MNRAS, 138, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Jørgensen, U. G., Hammer, D., Borysow, A., & Falkesgaard, J. 2000, A&A, 361, 283 [Google Scholar]
 Khomenko, E., Vitas, N., Collados, M., & de Vicente, A. 2017, A&A, 604, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Khomenko, E., Vitas, N., Collados, M., & de Vicente, A. 2018, A&A, 618, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kurucz, R. L. 1970, SAO Special Report, 309 [Google Scholar]
 Kurucz, R. L. 1993a, Atomic data for opacity calculations, CDROM No. 1 [Google Scholar]
 Kurucz, R. L. 1993b, Diatomic Molecular Data for Opacity Calculations, CDROM No. 15. [Google Scholar]
 Kurucz, R. L. 1993c, Opacities for Stellar Atmospheres: [+0.0],[+0.5],[+1.0], CDROM No. 2. [Google Scholar]
 Kurucz, R. L. 2002, AIP Conf. Ser., 636, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. L. 2014, in Determination of Atmospheric Parameters of B, A, F and GType Stars, eds. E. Niemczura, B. Smalley, & W. Pych (Berlin: Springer), 39 [CrossRef] [Google Scholar]
 Kurucz, R. L. 2017, Can. J. Phys., 95, 825 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. L., van Dishoeck, E. F., & Tarafdar, S. P. 1987, ApJ, 322, 992 [NASA ADS] [CrossRef] [Google Scholar]
 Labs, D. 1951, ZAp, 29, 199 [NASA ADS] [Google Scholar]
 Ludwig, H. G. 1992, PhD thesis, University of Kiel, Germany [Google Scholar]
 Ludwig, H. G., & Steffen, M. 2013, Mem. Soc. Astron. Ital. Suppl., 24, 53 [Google Scholar]
 Ludwig, H. G., Jordan, S., & Steffen, M. 1994, A&A, 284, 105 [NASA ADS] [Google Scholar]
 Magic, Z. 2014, PhD thesis, LudwigMaximilians University of Munich, Germany [Google Scholar]
 Magic, Z., Collet, R., Asplund, M., et al. 2013, A&A, 557, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nahar, S. N. 1995, A&A, 293, 967 [NASA ADS] [Google Scholar]
 Nordlund, A. 1982, A&A, 107, 1 [Google Scholar]
 Olson, G. L., & Kunasz, P. B. 1987, J. Quant. Spectrosc. Radiat. Transf., 38, 325 [CrossRef] [Google Scholar]
 Osorio, Y., Allende Prieto, C., Hubeny, I., Mészáros, S., & Shetrone, M. 2020, A&A, 637, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Palacios, A., Gebran, M., Josselin, E., et al. 2010, A&A, 516, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Peach, G. 1970, Mem. R. Astron. Soc., 73, 1 [Google Scholar]
 Pereira, T. M. D., Asplund, M., Collet, R., et al. 2013, A&A, 554, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rutten, R. J. 2003, Radiative Transfer in Stellar Atmospheres (Utrecht: Utrecht University Lecture Notes) [Google Scholar]
 Sbordone, L., Bonifacio, P., Castelli, F., & Kurucz, R. L. 2004, Mem. Soc. Astron. Ital. Suppl., 5, 93 [Google Scholar]
 Skartlien, R. 2000, ApJ, 536, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Strom, S. E., & Kurucz, R. L. 1966, J. Quant. Spectrosc. Radiat. Transf., 6, 591 [NASA ADS] [CrossRef] [Google Scholar]
 Trampedach, R., Asplund, M., Collet, R., Nordlund, Å., & Stein, R. F. 2013, ApJ, 769, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Vidal, C. R., Cooper, J., & Smith, E. W. 1973, ApJS, 25, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Vögler, A. 2003, PhD thesis, University of Göttingen, Germany [Google Scholar]
 Vögler, A., Bruls, J. H. M. J., & Schüssler, M. 2004, A&A, 421, 741 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Witzke, V., Shapiro, A. I., Cernetic, M., et al. 2021, A&A, 653, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wray, A. A., Ripoll, J. F., & Prabhu, D. 2006, in Studying Turbulence Using Numerical Simulation Database  XI., eds. P. Moin, & N. N. Mansour, Proceedings of the 2006 Summer Program, Center for Turbulence Research, Stanford University [Google Scholar]
These can be found in http://kurucz.harvard.edu/opacities.html. The particular table from ATLAS referred to in the present paper is in the section ‘Old Kurucz ODF files’, computed for the solar abundances from Anders & Grevesse (1989) and a turbulent velocity of 2 km s^{−1}. The corresponding monochromatic opacities are not available.
All Tables
Percentage of the combinations of the three separators (4 bins) with certain χ_{C} values for the four stellar types.
Percentage of the combinations of the three separators (4 bins) with certain χ_{H} values for the four stellar types.
Percentage of the combinations of the three separators (4 bins) with certain simultaneous χ_{C} and χ_{H} values for the four stellar types.
Deviation measures χ_{C} and χ_{H} for different binning strategies described in Sects. 5.1–5.3, for the four stars F3V, G2V, K0V, and M2V (columns are named after the corresponding sections).
All Figures
Fig. 1 Temperature stratification of the F3V (green), G2V (yellow). K0V (orange), and M2V (dark red) 1D hydrostatic models used in this work. The same colour coding is used throughout the paper to label results for the four stellar types. The grid used to solve the RTE is marked by dots. 

In the text 
Fig. 2 Monochromatic opacity (values shown in the colour–bars) for the F3V (first row), G2V (second row), KOV (third row), and M2V stars (fourth row) versus wavelength (horizontal axis) and geometrical height (vertical axis). Horizontal grey lines: Rosseland optical depths. The zero geometrical height is taken where logτ_{Ros} = 0. Other lines are shifted by Alogτ_{Ros} = 1. Black dashed thin line: height at which the continuum optical depth at each wavelength is unity (${\tau}_{\lambda}^{\text{cont}}=1$). Black solid thick line: heights at which the total optical depth τ_{λ} is unity. Values are averaged over 25 nm intervals. Tick marks on the top indicate wavelengths of the hydrogen series. Hydrogen lines gradually vanish with decreasing effective temperature. The colourbar values of each star range from min log ϰ to max log ϰ. 

In the text 
Fig. 3 Monochromatic radiative energy exchange rate Q_{v} (values shown in the colourbars) computed using the opacities from Fig. 2 for the models of the same four stellar types versus wavelength (horizontal axis) and geometrical height (vertical axis). Labels and overploted lines are same as in Fig. 2. The colourbar values of each star range from – max Q_{v} to max Q_{v}. Blue tones in the colourmap correspond to negative Q_{v}, i.e. cooling in the atmosphere, and red tones correspond to positive Q_{v}, i.e. heating. 

In the text 
Fig. 4 Three steps in construction of the ODF. Left: schematic representation of monochromatic opacities in an interval of wavelengths divided into steps of λ_{i}. The red vertical lines mark one ODF step with λ ∈ [λ_{2}, λ_{3}]. Middle: opacities of that step sorted by magnitude. The distribution is mapped onto the unit interval w ∈ [0, 1] and the wavelength dependence is lost within the step. Right: sorted opacities discretised in ODF substeps. 

In the text 
Fig. 5 Oscillator strength log gf versus wavelength for all molecular spectral lines included in SYNSPEC. Notice that the transitions that mainly contribute to the opacity have typically values of log gf > −5. 

In the text 
Fig. 6 Example of the radiative energy exchange rate Q versus geometrical height z for the M2V star. Top panel: the Q rate computed using the ODF from SYNSPEC data (black solid line) and the ODF from ATLAS (red solid line). The logarithm of the Rosseland optical depth is shown in the top axis. The area of the difference is shaded in grey. The zero of the geometrical height z is at the continuum optical depth at 500 nm τ_{5} = 1. The range of heights of the cooling (Q < 0) and heating (Q > 0) components are indicated in cyan and orange respectively. The heights z_{b}, z_{c→h}, and z_{t} are used as integral bounds to calculate the deviation measures χ_{C} and χ_{Η} (Eqs. (13) and (14)). The height z_{c→h} is also marked with a vertical black line. Bottom panel: absolute difference between the Q values. 

In the text 
Fig. 7 Bolometric radiative energy exchange rate Q^{ODF} against the geometrical height in the four stellar models. The results using ODF from SYNSPEC (coloured lines) are compared to those using the ODF from ATLAS (thin black line). Top left of each panel: scaling factor of vertical axis. Bottom left of each panel: difference between the two ODFs computed using Eqs. (13) and (14). Top right of each panel: ratio of the heating area with respect to the cooling area for the Q for both data sets (A for ATLAS and S for SYNSPEC). In cyan in the bottom right panel for the M2V star: ODF computed using monochromatic opacity from SYNSPEC, excluding the lines from TiO and VO. The corresponding deviations for that Q curve are χ_{C} ≃ 7 and χ_{Η} ≃ 29. 

In the text 
Fig. 8 Example of the reference optical depth and wavelength of the ODF points organised into bins for the F3V, G2V, K0V, and M2V stars. As explained in the text, the Rosseland optical depth shown in the vertical axis is used as the reference τ^{ref}(z_{i,j}), where z_{i,j} is the geometrical height at which τ_{i,j} = 1 for each step and substep of the ODF. We show the same {τ, λ} binning configuration for all the stars. The colours indicate different bins. Each point corresponds to one ODF substep, all centred in the middle wavelength of their corresponding step. The last column of dots at the right of each panel corresponds to the ODF step with λ ∈ [9600, 9.5 × 10^{4}] nm, and thus, appears separated from the rest of the dots. 

In the text 
Fig. 9 Histogram of the logarithm of the relative number of points from the ODF as a function of τ^{ref}(z_{i,j}) in the vertical axis and λ in the horizontal, as in Fig. 8. The density of the ODF points is the largest in the continuum (especially in the visible) and in the UV owing to the sampling of the ODF steps (see the description of the ODF tables in Sect. 4). 

In the text 
Fig. 10 Deviation measures χ_{C} (upper row) and χ_{Η} (lower) for the G2V star as function of the location of the separator τ_{1} when the other two separators are fixed. Columns correspond to six different fixed values of τ_{0} (values indicated at the top of each column). Each curve corresponds to a different location of the top separator τ_{2}. The colour of each curve indicates the location of τ_{2} used to compute that curve. We do not show all the possible locations for τ_{0} and τ_{2} to avoid overcrowding of the plots. The grey horizontal lines show χ_{C} = 4% and χ_{Η} = 15%. 

In the text 
Fig. 11 Same as in Fig. 10, but for the M2V star. The grey horizontal lines show χ_{C} = 9% and χ_{H} = 27%. 

In the text 
Fig. 12 Deviation measures χ_{C} (upper panel) and χ_{H} (lower) versus the number of bins for the four stars (colours as in Fig. 1). Solid lines: deviation measures for the Q computed using 3, 4, , 11, 21, 31, 41, 51, and 101 τbins distributed equidistantly. Circles: deviations for Q computed with the Rosseland opacity. Pluses: the deviation of the optimal combinations of four τbins (see Sect. 5.1 and Table 4). 

In the text 
Fig. 13 Example of the binning used to compute the Q values for the G2V star, using different numbers of {τ, λ}bins. Similarly to Fig. 8, each grey dot corresponds to one ODF substep. For this star, the τseparators in the UV in all the experiments are fixed at log τ = {−1.8, −4.2}. Columns, from left to right, show the cases for N = 6–9 bins. The red numbers show indices of the bins. The case for N = 6 has two λseparators at 400 and 1600 nm. For more than 6 bins, bins 1 to 5 are kept identical, but more λ separators are added between 400 and 9750 nm. A similar binning setup is used for the other stars with modified τseparators in the UV and the λ separator between the short and long wavelengths for N > 6. 

In the text 
Fig. 14 Deviation measures χ_{C} (upper panel) and χ_{Η} (lower) versus the number of bins, for the four stars (colours as in Fig. 1), for the case of {τ, λ)binning approach similar to Magic (2014). Solid lines: deviation measures for the Q computed using 6, 7, 8, 9, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, and 100 bins distributed as in Fig. 13. Circles: deviations for Q computed with the Rosseland opacity. Crosses: deviations for Q computed with only two τseparators located at the heights of the bins in the UV (no λseparators). Pluses: the deviation of the optimal combinations of four τbins (see Sect. 5.1 and Table 4). 

In the text 
Fig. 15 Comparison of the radiative energy exchange rate Q from bin 7 for N = 7 and bins 8 and 9 for N = 9 corresponding to the cases shown in Fig. 14 for the M2V star. Left panel: Q_{7} dominates the total Q^{OB} = Σ_{l} Q_{l} for N = 7. Right panel: comparison of Q_{8} and Q_{9} for N = 9 against Q_{7} for N = 7. 

In the text 
Fig. 16 Deviation measures χ_{C} (upper panel) and χ_{Η} (lower) against the number of bins for the four stars, computed for the strategy with five fixed λseparators at λ ∈ (381, 726, 1383, 2636, 5023) nm and with varying number of τseparators (1, 2, 3, 5, 6, 8, 10, 12, 14, and 16) and therefore varying total number of bins (12, 18, 24, 36, 42, 54, 66, 78, 90, and 102 bins). Solid lines: deviation measures for Q computed using a different number of {τ, λ}bins. Circles: deviations for Q computed with the Rosseland opacity. Pluses: deviations of the optimal combinations of four τbins (see Sect. 5.1 and Table 4). 

In the text 
Fig. 17 Comparison of Q and Q/ρ for the four stars (each row) and five opacity approaches: ODF (black solid line), grey opacity computed with Eq. (C.1; grey solid line), and the binned opacities for four bins (red solid line), 15 bins (yellow solid line), and 24 bins (black dotted line) corresponding to Cols. 5.1, 5.3a and 5.3b in Table 4. Left column: Q for the deeper part of the atmosphere where mainly cooling occurs. Second column: Q/ρ for the same part of the atmosphere. Third column: Q for the higher layers where the heating occurs. Right column: Q/ρ for the same part of the atmosphere. 

In the text 
Fig. B.1 ODF opacity against the temperature of the four stellar atmospheres (different colour lines, corresponding to a different spectral type as in Fig. 1), and for 12 substeps (different panels) of the ODF step defined by the λ ∈ [590, 600] nm range: SYNPEC (thick coloured lines) and ATLAS (thin black lines). The number at the left bottom corner of each panel indicates the substep in the order of increasing opacities. This ODF step is an example of a generally good match between SYNPEC and ATLAS that is found in the majority of ODF steps in the visible and IR. The only systematic difference is present in all substeps at the top part of the M model. 

In the text 
Fig. B.2 Same as Fig. B.1, but for the step defined by λ ∈ [650, 660] nm containing hydrogen Ha. This example shows that there are significant differences between the two data sets in the strong hydrogen line. The differences are present in all substeps for the M star, but also in the more opaque substeps (rows three and four) for the three hotter stellar types. 

In the text 
Fig. B.3 Same as Fig. B.1, but for the step defined by λ ∈ [580, 590] nm containing NaI D1 and D2 line. 

In the text 
Fig. B.4 Same as Fig. B.1, but for the step defined by λ ∈ [193, 197.5] nm. The differences in the UV between the two tables are large in all the substeps and for all stellar models. 

In the text 
Fig. B.5 Radiative energy exchange rate Q_{i} computed from the four model atmospheres using the opacity shown in Fig. B.1 and integrated over the interval λ ∈ [590, 600] nm. The values on the vertical axis are scaled with the factor indicated in the top left of each panel. Bottom left of each panel: deviation between the two ODFs computed using Eqs. 13, 14. Top right of each panel: ratio of the heating area with respect to the cooling area for the Q_{i} for both data sets (A for ATLAS and S for SYNSPEC). 

In the text 
Fig. B.6 Same as Fig. B.5, but for the λ ∈ [650, 660] nm interval that contains the hydrogen Hα line. Computed using the opacity shown in Fig. B.2. 

In the text 
Fig. B.7 Same as Fig. B.5, but for λ ∈ [580, 590] nm, interval that contains the NaI D1 and D2 line. Computed using the opacity shown in Fig. B.3. 

In the text 
Fig. B.8 Same as Fig. B.5, but for λ ∈ [193, 197.5] nm, computed using the opacity shown in Fig. B.4. 

In the text 
Fig. C.1 Distribution of the opacity bins in optical depth and wavelength for the four stars. Similarly to Fig. 8, each coloured dot corresponds to one ODF substep. The wavelength separators are identical for all stars (λ < 300, ∈ [300, 500], ∈ [500, 700], ∈ [700, 1000], ∈ [1000 – 2000], and > 2000 nm). All values within one wavelength range are shown in the same colour. The τseparators are selected from the optimal combinations found in Sect. 5.1 for each star: log τ = {−1.9, −5, −7} for the F3V star; log τ = {0, −2.5, −5} for the G2V; log τ = {−0.2, −3, −5} for the K0V; log τ = {0.3, −0.8, −2.5} for the M2V. 

In the text 
Fig. C.2 Opacity ${\overline{\varkappa}}_{l}$ of each τbin in pairs of λranges for the four stars (spectral type at the top left of each panel) versus the geometrical height (bottom axis of each panel) and optical depth (top axis). The λranges are selected to illustrate the different variations of the opacity with height depending on the spectral region. The colours of the curves correspond to the λranges in Fig. C.1, with the values of each range in the legend. 

In the text 
Fig. C.3 Opacity ${\overline{\varkappa}}_{l}$ for each τbin in λ < 300 nm computed with Eq. 20 (solid blue curves) compared to the opacity ${\overline{\varkappa}}_{l}={\overline{\varkappa}}_{l}^{\text{Ro}}$ (dotted black curves) for the F3V star. The rest of the coloured curves correspond to the opacity for the least opaque bin for each wavelength range. Dashed black line: Rosseland opacity mean computed with Eq. C.1. The bottom and top axis show the geometrical height and optical depth, respectively. 

In the text 
Fig. C.4 Rosseland opacity mean computed with Eq. C.1 (black dashed line) and Eq. C.2 (red thin solid line) for the F3V star. The bottom and top axis show the geometrical height and optical depth, respectively. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.