11institutetext: Instituto de Astrofísica de Canarias, 38200 La Laguna,
Tenerife, Spain
11email<EMAIL_ADDRESS>22institutetext: Departamento de Astrofísica de
la Universidad de La Laguna, 38200 La Laguna, Tenerife, Spain 33institutetext:
Isaac Newton Group of Telescopes, Apartado de Correos 321, 38700 Santa Cruz de
La Palma, Canary Islands, Spain 44institutetext: Steward Observatory,
University of Arizona, Tucson, USA

# Opacity for realistic 3D MHD simulations of cool stellar atmospheres

A. Perdomo García , 1122 N. Vitas, 1122 E. Khomenko , 1122 M. Collados , 1122
C. Allende Prieto, 1122 I. Hubeny 44 Y. Osorio ,, 112233

(Received Month day, year; accepted Month day, year)

###### Abstract

Context. Realistic 3D time-dependent simulations of stellar near-surface
convection employ the opacity binning method for efficient and accurate
computation of the radiative energy exchange. The method provides several
orders of magnitude of speed-up, 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 are used
to construct opacity distribution function (ODF) that is then verified through
detailed comparison with the results of the ATLAS code. The opacity binning
method is implemented with the SYNSPEC opacities for four representative cool
main-sequence 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 case of the F, G,
and K – type stars. Significant differences, coming mainly from the molecular
line lists, are found for the M – type 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 K – type stars. In
the case of the M – type 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

##  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: the absorption coefficients of all contributing
(bound-bound, bound-free, and free-free) 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 late-type 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 $\approx$$10^{9}$ lines and
$\approx$$10^{6}$ wavelength points, is difficult even in one-dimensional (1D)
modelling, and virtually impossible in three-dimensional (3D) time-dependent
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 $\rho$); 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 NLTE lines, which are computed explicitly. The opacities are computed as
by-products 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 wavelength-integrated source term $Q$. The distribution of radiative
heating and cooling with wavelength is irrelevant, at least as long as the
one-fluid 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 monotonic distribution that is then
discretised by averaging over several ($\approx$$10$) substeps. In that way
the RTE can be solved for several thousands ODF values instead 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 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 time-dependent 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 main-sequence 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 up-to-date 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
source111https://www.as.arizona.edu/ehubeny/tlusty208-package/tl208-s54.tar.gz
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
near-surface convection in cool stars (Perdomo et al, in preparation) 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 1D 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.

##  2 Method

###  2.1 Model atmospheres

Figure 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 4 stellar
types. The grid used to solve the RTE is marked by dots.

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(\tau)$ (see
chapter 17 from Hubeny & Mihalas 2014, and references 507–510 there) of the
the Harvard-Smithsonian Reference Atmosphere (Gingerich et al. 1971) with the
effective temperature of each star, following the expression
$T^{4}(\tau)=\frac{3}{4}T_{\mathrm{eff}}^{4}\times\left[\tau+q(\tau)\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 plane-parallel static atmosphere in LTE without considering line
scattering is

$\mu\frac{dI_{\nu}(z,\nu)}{dz}=\rho(z)\varkappa_{\nu}(z)\left[B_{\nu}(z)-I_{\nu}(z)\right],$
(1)

where $I_{\nu}(z,\nu)$ is the specific intensity (erg s-1 cm-2 Hz-1 ster-1) at
frequency $\nu$ in the direction $\mu$ and at the geometrical height $z$ in
the atmosphere; $B_{\nu}$ is the Planck function; $\rho(z)$ is the mass
density; and $\varkappa_{\nu}(z)$ is the monochromatic absorption coefficient
per unit mass (cm2 g-1; referred to as opacity throughout the text). The
opacity is computed as the sum of the continuum opacity (mainly contributions
from the free-free and bound-free processes) and the line opacity (bound-bound
individual transitions included in the line lists for atoms and molecules),

$\varkappa_{\nu}=\varkappa_{\nu}^{\mathrm{cont}}+\varkappa_{\nu}^{\mathrm{lines}}$
(2)

We restrict ourselves to the LTE case that is valid only in deep layers where
$J\approx B$. Therefore our treatment of the coherent scattering terms (see
App. A) in the total opacity is only approximate (e.g. Hayek 2010).

We solve the RTE in the vertical direction ($\mu=\pm 1$). At each frequency
$\nu$, the formal solution for outward and inward intensity ($I^{+}$, $I^{-}$)
in LTE is:

$\begin{split}I^{\pm}_{\nu}(z_{i})&=I^{\pm}_{\nu}(z_{i\mp
1})\exp\left(-\Delta\tau_{\nu}^{i\mp
1}\right)+\int_{0}^{\Delta\tau_{\nu}^{i\mp
1}}B_{\nu}(t)\exp\left(t-\Delta\tau_{\nu}^{i\mp 1}\right)dt\equiv\\\ &\equiv
I^{\pm}_{\nu}(z_{i\mp 1})\exp\left(-\Delta\tau_{\nu}^{i\mp 1}\right)+\Delta
I^{\pm}_{\nu}(z_{i}),\end{split}$ (3)

where $i$ is the height index for each point in the atmosphere, being zero at
the bottom and increasing upwards; $\tau_{\nu}^{i}$ is the optical depth at
the height with index $i$ ($d\tau_{\nu}(z)=\rho(z)\varkappa_{\nu}(z)dz$); and
$\Delta\tau_{\nu}^{i+1}=\tau_{\nu}^{i+1}-\tau_{\nu}^{i}$ and
$\Delta\tau_{\nu}^{i-1}=\tau_{\nu}^{i}-\tau_{\nu}^{i-1}$. To solve the RTE, we
apply the short characteristics method (Olson & Kunasz 1987), using the linear
approximation of the Planck function:

$\Delta I_{\nu}^{\pm}(z_{i})=\psi_{0}B_{\nu}(z_{i})+\psi_{1}B_{\nu}(z_{i\mp
1}).$ (4)

For each frequency $\nu$ we compute the coefficients for the local and the
upwind points

$\begin{split}\psi_{0}&=1-\frac{1}{\Delta\tau_{\nu}}+\frac{1}{\Delta\tau_{\nu}}\exp\left(-\Delta\tau_{\nu}\right),\\\
\psi_{1}&=\frac{1}{\Delta\tau_{\nu}}-\frac{\Delta\tau_{\nu}+1}{\Delta\tau_{\nu}}\exp\left(-\Delta\tau_{\nu}\right).\\\
\end{split}$ (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 $\Delta\tau_{\nu}$ in
each case ($\Delta\tau_{\nu}=\Delta\tau_{\nu}^{i-1}$ for the upward
intensities; $\Delta\tau_{\nu}=\Delta\tau_{\nu}^{i+1}$ for the inward ones).
This first order short-characteristics 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_{\nu}$ (erg s-1 cm-2 Hz-1 ster-1) and the flux
$\@vec{\mathcal{F}}_{\nu}$ (erg s-1 cm-2 Hz-1) can then be calculated for a
vertical ray from the known $I_{\nu}$:

$J_{\nu}(z_{i})=\frac{1}{2}\int_{-1}^{1}I_{\nu}(z_{i})d\mu\equiv\frac{1}{2}\left(I^{+}_{\nu}(z_{i})+I^{-}_{\nu}(z_{i})\right),$
(6) $\mathcal{F}_{\nu}(z_{i})=2\pi\int_{-1}^{1}I_{\nu}(z_{i})\mu d\mu\equiv
2\pi\left(I^{+}_{\nu}(z_{i})-I^{-}_{\nu}(z_{i})\right).$ (7)

Finally, the monochromatic radiative energy exchange rate $Q_{\nu}$ (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^{F}_{\nu}(z_{i})=-\@vec{\nabla}\@vec{\mathcal{F}}_{\nu}\equiv-(\mathcal{F}_{\nu}(z_{i+1})-\mathcal{F}_{\nu}(z_{i-1}))/(z_{i+1}-z_{i-1}),$
(8)

or directly from the mean intensity,

$Q^{J}_{\nu}(z_{i})=4\pi\varkappa_{\nu}(z_{i})\rho(z_{i})\left(J_{\nu}-B_{\nu}\right).$
(9)

The monochromatic rate is calculated as a weighted combination of the two,
using equation 4.13 from Vögler (2003) (following the suggestion from Bruls et
al. 1999):

$Q_{\nu}(z_{i})=\exp\left(-\tau^{i}_{\nu}/\tau_{\mathrm{h}}\right)Q^{J}_{\nu}(z_{i})+\left[1-\exp\left(-\tau^{i}_{\nu}/\tau_{\mathrm{h}}\right)\right]Q^{F}_{\nu}(z_{i}),$
(10)

where $\tau_{\mathrm{h}}=0.1$. As it is explained in Vögler (2003), this
smooth transition between $Q^{J}_{\nu}$ and $Q^{F}_{\nu}$ avoids the accuracy
problem that $Q^{J}_{\nu}$ has in the optically thick regime (where $J_{\nu}$
approaches $B_{\nu}$) 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^{F}_{\nu}$.

##  3 Monochromatic opacity

Figure 2: Monochromatic opacity (values shown in the colour–bars) for the F3V
(1st row), G2V (2nd row), K0V (3rd row), and M2V stars (4th row) versus
wavelength (horizontal axis) and geometrical height (vertical axis).
Horizontal grey lines: Rosseland optical depths. The zero geometrical height
is taken where $\log\tau_{\mathrm{Ros}}=0$. Other lines are shifted by
$\Delta$log$\tau_{\mathrm{Ros}}=1$. Black dashed thin line: height at which
the continuum optical depth at each wavelength is unity
($\tau^{\mathrm{cont}}_{\lambda}=1$). Black solid thick line: heights at which
the total optical depth $\tau_{\lambda}$ 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 colour-bar values of each star range from $\min\log\varkappa$ to
$\max\log\varkappa$. Figure 3: Monochromatic radiative energy exchange rate
$Q_{\nu}$ (values shown in the colour–bars) 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 colour-bar values of each star
range from $-\max|Q_{\nu}|$ to $\max|Q_{\nu}|$. Blue tones in the colour-map
correspond to negative $Q_{\nu}$, i.e. cooling in the atmosphere, and red
tones correspond to positive $Q_{\nu}$, i.e. heating.

For the computation of the monochromatic opacity $\varkappa_{\nu}$ 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 the 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 preparation), 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, and Modestov et al. in preparation).

The EOS in SYNSPEC is solved for a specified grid of temperature $T$ and mass
density $\rho$. 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 (2017a, b).

Once the EOS is solved, the opacity is computed in the same ($T,\rho$) grid,
using the atomic and molecular data from several sources (see App. A). The
data used by the code can be customised by the user by changing, for example,
the line lists for the bound-bound 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 App. 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 we will be comparing to. The $T$-axis of the grid has 60 values
in the range $T\in[1995,125900]$ K, with constant step $\Delta\log T\simeq
0.03$ K (in the present text, $\log\equiv\log_{10}$). The $\rho$-axis has 35
values in the range $\rho\in[10^{-13},3.2\times 10^{-3}]$ g cm-3 with
$\Delta\log\rho\simeq 0.3$. The table is computed for wavelengths in the range
$\lambda\in[20,9.5\times 10^{4}]$ nm, with $\Delta\log\lambda=10^{-6}$
($\lambda$ in nm). In total there is around $7.7\times 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, $T\in[1995,19900]$ K and
$\rho\in[10^{-9},3.2\times 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$ nm
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 $\tau_{\nu}=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 sub-intervals 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^{\mathrm{cont}}_{\nu}=1$. The ionisation edges and the typical shape of
the bound-bound and bound-free H- opacity coefficient are visible for the
$\tau^{\mathrm{cont}}_{\nu}=1$ curve. In the UV the line forest elevates the
height of the $\tau_{\lambda}=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
$\tau_{\lambda}=1$ curve is shaped by millions of molecular transitions.

Following the method from Sect. 2, we solve the RTE to obtain $Q_{\nu}$. Once
$Q_{\nu}$ is computed, we can integrate it in frequency using the trapezoidal
rule:

$Q=\int_{\nu}Q_{\nu}d\nu\equiv\frac{1}{2}\sum_{i}(Q_{\nu_{i}}+Q_{\nu_{i+1}})\Delta\nu_{i}.$
(11)

The computed $Q_{\nu}$ 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_{\nu}$ is large for a wider wavelength range and the significance
of the contribution to $Q_{\nu}$ 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_{\nu}$ for the heating in the atmosphere
($Q_{\nu}>0$) and the cooling ($Q_{\nu}<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 near-
UV and blue wavelengths. For the M star, the heating is comparable to the
cooling in magnitude, and ranges from the visible to the near-IR. 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/\rho$ (erg g-1 s-1), which would put more
emphasis on the radiation exchange at lower densities.

##  4 Opacity distribution functions

###  4.1 Construction of the opacity distribution function

Figure 4: Three steps in construction of the ODF. Left: schematic
representation of monochromatic opacities in an interval of wavelengths
divided into steps of $\lambda_{i}$. The red vertical lines mark one ODF step
with $\lambda\in[\lambda_{2},\lambda_{3}]$. Middle: opacities of that step
sorted by magnitude. The distribution is mapped onto the unit interval
$w\in[0,1]$ and the wavelength dependence is lost within the step. Right:
sorted opacities discretised in ODF substeps. Figure 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$. Figure 6: Top panel: The radiative
energy exchange rate $Q$ versus geometrical height $z$ for the M2V star
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
$\tau_{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_{\mathrm{b}}$, $z_{\mathrm{c\rightarrow h}}$, and $z_{\mathrm{t}}$ are used
as integral bounds to calculate the deviation measures $\chi_{\mathrm{C}}$ and
$\chi_{\mathrm{H}}$ (Eqs. 13, 14). The height $z_{\mathrm{c\rightarrow h}}$ is
also marked with a vertical black line. Bottom panel: absolute difference
between the $Q$ values. Figure 7: Bolometric radiative energy exchange rate
$Q^{\mathrm{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, 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
$\chi_{\mathrm{C}}\simeq 7$ and $\chi_{\mathrm{H}}\simeq 29$.

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 time-
dependent 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
$\lambda\in[\lambda_{2},\lambda_{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 $\omega_{j}$
($\sum_{j}\omega_{j}=1$). The length of a substep $\Delta\lambda_{i,j}$ is
calculated in terms of the length of the step $\Delta\lambda_{i}$ and the
weight $\omega_{j}$, so that
$\Delta\lambda_{i,j}=\omega_{j}\Delta\lambda_{i}$222The length
$\Delta\lambda_{i,j}$ is not an interval in the wavelength scale, but it
counts a certain number of points in the distribution.. 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 $\rho$ (as in
SYNSPEC333SYNSPEC can also work with a grid of electron density or pressure.)
or the gas pressure $P$ (as in the ATLAS code, Kurucz 1970). Kurucz (1993c)
created ODF tables444These 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. 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,\rho$) 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, $\lambda\in[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 $\lambda$ 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
$\omega_{j}=\left\\{0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.05,2/60,1/60\right\\}$.
This kind of non-uniform 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 cm2 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 (H2, HD, MgH, NH, CH, SiH, OH, CH, C2, 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 App. A and Fig. 5) and includes an up-to-date 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.

###  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 $\rho$ along the four 1D model of atmospheres. Four examples
representing typical differences are shown in the App. B, in Figs. 18-21
(where black lines represent the opacity from ATLAS along the models, and the
coloured lines on top of them represent the 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^{\mathrm{ODF}}=\sum_{i}Q_{i}\equiv\sum_{i}\Delta\lambda_{i}\sum_{j}Q_{i,j}\omega_{j},$
(12)

where $Q_{i}$ is the rate integrated in wavelength in the $i$th 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 $\tau=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:

$\mathrm{\chi_{\mathrm{C}}}=\frac{A\left(\left|Q^{(1)}-Q^{(2)}\right|\right)_{\mathrm{C}}}{A\left(\left|Q^{(1)}\right|\right)_{\mathrm{C}}},$
(13)
$\mathrm{\chi_{\mathrm{H}}}=\frac{A\left(\left|Q^{(1)}-Q^{(2)}\right|\right)_{\mathrm{H}}}{A\left(\left|Q^{(1)}\right|\right)_{\mathrm{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(z)\right)_{C}=\int_{z_{\mathrm{b}}}^{z_{\mathrm{C\rightarrow
H}}}f(z)dz$ and $A\left(f(z)\right)_{H}=\int_{z_{\mathrm{C\rightarrow
H}}}^{z_{\mathrm{t}}}f(z)dz$. The height at $z_{\mathrm{C\rightarrow H}}$ is
the height where $Q$ changes sign for the first time above the cooling
component; $z_{\mathrm{b}}$ and $z_{\mathrm{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 $\left|Q\right|<2\times
10^{-4}\left|\min(Q)\right|$ (see Fig. 6). If
$A\left(Q\right)_{\mathrm{H}}<1\%\;A\left(Q\right)_{\mathrm{C}}$,
$\chi_{\mathrm{H}}$ is not computed (these cases are marked as
$\chi_{\mathrm{H}}=$ \- - - $\%$ in the figures).

The radiative energy exchange terms $Q_{i}$ computed from the four models
using the same steps as in Figs. 18-21 are shown in App. B, in Figs. 22-25.
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
($\lambda\in[590,600]$ nm) is shown in Fig. 18. The $Q_{i}$ values computed in
the same step closely overlap for these stars with all $\chi_{\mathrm{C}}$ and
$\chi_{\mathrm{H}}$ between $1.4\%$ and $3.1\%$ (see, for example, Fig. 22).
For the near-UV, visible, and IR, the steps that include strong lines show
significant differences in the opacity. Some examples are the steps that
include hydrogen H$\alpha$ (Fig. 19) or NaI D1 and D2 (Fig. 20). Similar
discrepancies are found in the steps containing other strong lines like the
CaII H and K lines, the G band, HI $\beta$, 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. 19). 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. 20). The
number of affected substeps varies between the 2 and 4 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 $\chi_{\mathrm{C}},\chi_{\mathrm{H}}\lesssim
5-10\%$ for the F, G, and K stars (see, for example, Figs. 23,24). An
exception happens for the G star for the step containing hydrogen H$\alpha$
where $\chi_{\mathrm{H}}\simeq 25\%$ (see Fig. 23). 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. 21). These
discrepancies are also evident in the corresponding values of $Q_{i}$ (see an
example in Fig. 25), for which we find differences of around
$\chi_{\mathrm{C}},\chi_{\mathrm{H}}\simeq 5-40\%$ for $\lambda\gtrsim 190$ nm
(as in the example from Fig. 25), and around
$\chi_{\mathrm{C}},\chi_{\mathrm{H}}\simeq 30-80\%$ for $\lambda<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 cross-sections from the Opacity Project
555http://cdsweb.u-strasbg.fr/topbase/publi.html 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 cross-sections 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
$\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ are also larger, for most of the
steps in the visible and NIR they are in the range $\simeq$10–35%
($\simeq$40–60% in the worst cases).

As previously discussed, the discrepancies between the two tables for the
later-type 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 $\chi_{\mathrm{H}}\simeq 77$%, in contrast
to the values $\chi_{\mathrm{C}}<3\%$ and $\chi_{\mathrm{H}}<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 $\chi_{\mathrm{H}}$
is considerably reduced, but it is still relatively high ($\simeq$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 H2O and TiO, from EXOMOL
666http://www.exomol.com/data, or different partition functions and oscillator
strengths.

##  5 Opacity binning

Figure 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 $\tau^{\mathrm{ref}}(z_{i,j})$, where $z_{i,j}$ is the
geometrical height at which $\tau_{i,j}=1$ for each step and substep of the
ODF. We show the same $\\{\tau,\lambda\\}$ 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 $\lambda\in[9600,9.5\times 10^{4}]$ nm, and thus, appears separated from
the rest of the dots. Figure 9: Histogram of the logarithm of the relative
number of points from the ODF as a function of $\tau^{\mathrm{ref}}(z_{i,j})$
in the vertical axis and $\lambda$ 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).

A further approximation of the opacity computation is the OB method initially
introduced by Nordlund (1982). Its implementation in the MURaM and CO5BOLD
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 $\tau^{\mathrm{ref}}$ of a representative one-
dimensional model atmosphere. We use the Rosseland optical depth as reference.
An opacity bin $l$ is defined by selecting a pair of lower and upper
$\tau$-separators ($\tau^{\mathrm{ref,l}}_{0}$, $\tau^{\mathrm{ref,l}}_{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 ($\tau$-separators) and one in wavelength
($\lambda$-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}(z)=\tau_{i,j}(z_{\mathrm{top}})\,+\int_{z_{\mathrm{top}}}^{z}\varkappa_{i,j}(z)\,\rho\,\mathrm{d}z,$
(15)

where $\rho$ 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
$\varkappa_{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^{\mathrm{ref,l}}_{0}\leq\tau^{\mathrm{ref}}(z_{i,j})<\tau^{\mathrm{ref,l}}_{1}$,
where $z_{i,j}$ is the height at which $\tau_{i,j}=1$. Additionally, when the
$\lambda$-separators are used, every ODF value is sorted by wavelength so that
the central wavelength of an ODF step falls between a pair of
$\lambda$-separators. Figure 8 shows examples of such grouping of the ODF
values into five bins for our four stellar models. There are two
$\lambda$-separators approximately dividing the wavelength axis into
ultraviolet, visible, and infrared and there is one $\tau$-separator for the
UV ($\log\tau=-3$) and a joint one for the visible and IR ($\log\tau=-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_{\mathrm{top}}$ to all ODF values with
$\tau_{i,j}(z_{\mathrm{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,\rho$ of the grid of the ODF we compute the Planck
function,

$\mathrm{B}_{l}=\sum_{i(l)}\Delta\lambda_{i}\mathrm{B}_{i},$ (16)

its derivative,

$\left.\frac{\partial\mathrm{B}}{\partial\mathrm{T}}\right|_{l}=\sum_{i(l)}\Delta\lambda_{i}\frac{\partial\mathrm{B}_{i}}{\partial\mathrm{T}},$
(17)

the Planck mean opacity,

$\overline{\varkappa}^{\mathrm{Pl}}_{l}=\left(\sum_{i(l)}\Delta\lambda_{i}\mathrm{B}_{i}\sum_{j(i,l)}\omega_{j}\varkappa_{i,j}\right)\Bigg{/}\mathrm{B}_{l},$
(18)

and the Rosseland mean opacity,

$\overline{\varkappa}^{\mathrm{Ro}}_{l}=\left.\frac{\partial\mathrm{B}}{\partial\mathrm{T}}\right|_{l}\Bigg{/}\left(\sum_{i(l)}\Delta\lambda_{i}\frac{\partial\mathrm{B}_{i}}{\partial\mathrm{T}}\sum_{j(i,l)}\frac{\omega_{j}}{\varkappa_{i,j}}\right),$
(19)

where $\Delta\lambda_{i}$, $\omega_{j}$, and $\varkappa_{i,j}$ are,
respectively, the length of the wavelength step, the weight of the substep,
and the opacity of the $i$th step and $j$th substep of the ODF.
$\mathrm{B}_{i}$ is the Plank function in the middle of the wavelength step
for each temperature $T$ in the ($T,\rho$) grid used to compute the ODF.

Finally, following Ludwig (1992), $\overline{\varkappa}^{\mathrm{Pl}}_{l}$ and
$\overline{\varkappa}^{\mathrm{Ro}}_{l}$ are combined to get a mean value of
the opacity for each bin:

$\overline{\varkappa}_{l}=\left(2^{-\tau_{l}/\tau_{\mathrm{thr}}}\right)\overline{\varkappa}^{\mathrm{Pl}}_{l}+\left(1-2^{-\tau_{l}/\tau_{\mathrm{thr}}}\right)\overline{\varkappa}^{\mathrm{Ro}}_{l},$
(20)

where the transition between thin and thick optical regimes is at the
threshold optical depth $\tau_{\mathrm{thr}}=0.35$ and $\tau_{l}$ is computed
using pressure $p$ in the ODF grid and a stellar surface gravity $g$,
approximately as $\tau_{l}=\overline{\varkappa}^{\mathrm{Ro}}_{l}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. 28 in App. C).

As in the ODF case, we compute the bolometric $Q$ rate, $Q^{\mathrm{OB}}$, by
solving the RTE for each bin and summing over the bins:

$Q^{\mathrm{OB}}=\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 $\tau$ and $\lambda$, 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 height-dependent absolute difference of
$Q^{(1)}$ and $Q^{(2)}$ divided by the density stratification $\rho$ 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 trial-and-error 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 free-parameters, while it preserves
separate information on both cooling and the heating component of $Q$. As this
is a multi-parameter study we break the problem down by first studying the
influence of the location of the $\tau$-separators for a fixed number of bins
(4) and with no bins in the wavelength direction. We then allow the number of
$\tau$ bins to change. Finally, we test several cases including the $\lambda$
bins.

###  5.1 $\tau$-bin location

Table 1: Percentage of the combinations of the three separators (4 bins) with certain $\chi_{\mathrm{C}}$ values for the four stellar types. | Percentage of combinations
---|---
| $\chi_{\mathrm{C}}<20\%$ | $\chi_{\mathrm{C}}<10\%$ | $\chi_{\mathrm{C}}<5\%$ | $\chi_{\mathrm{C}}<3\%$
F | 98 | 74 | 2.3 | 0
G | 100 | 96 | 87 | 43
K | 100 | 99 | 93 | 34
M | 100 | 28 | 2.4 | 0
Table 2: Percentage of the combinations of the three separators (4 bins) with certain $\chi_{\mathrm{H}}$ values for the four stellar types. | Percentage of combinations
---|---
| $\chi_{\mathrm{H}}<50\%$ | $\chi_{\mathrm{H}}<30\%$ | $\chi_{\mathrm{H}}<20\%$ | $\chi_{\mathrm{H}}<10\%$
F | 64 | 55 | 50 | 16
G | 85 | 78 | 55 | 8.9
K | 83 | 77 | 62 | 6.8
M | 56 | 7.4 | 0.4 | 0
Table 3: Percentage of the combinations of the three separators (4 bins) with certain simultaneous $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ values for the four stellar types. | Percentage of combinations
---|---
| $\chi_{\mathrm{C}}<10\%$ | $\chi_{\mathrm{C}}<10\%$ | $\chi_{\mathrm{C}}<5\%$ | $\chi_{\mathrm{C}}<5\%$
| $\chi_{\mathrm{H}}<50\%$ | $\chi_{\mathrm{H}}<20\%$ | $\chi_{\mathrm{H}}<20\%$ | $\chi_{\mathrm{H}}<15\%$
F | 64 | 50 | 1.8 | 1.2
G | 85 | 55 | 55 | 32
K | 83 | 62 | 62 | 25
M | 14 | 0.4 | 0 | 0
Figure 10: Deviation measures $\chi_{\mathrm{C}}$ (upper row) and
$\chi_{\mathrm{H}}$ (lower) for the G2V star as function of the location of
the separator $\tau_{1}$ when the other two separators are fixed. Columns
correspond to six different fixed values of $\tau_{0}$ (values indicated at
the top of each column). Each curve corresponds to a different location of the
top separator $\tau_{2}$. The colour of each curve indicates the location of
$\tau_{2}$ used to compute that curve. We do not show all the possible
locations for $\tau_{0}$ and $\tau_{2}$ to avoid overcrowding of the plots.
The grey horizontal lines show $\chi_{\mathrm{C}}=4\%$ and
$\chi_{\mathrm{H}}=15\%$. Figure 11: Same as in Fig. 10, but for the M2V star.
The grey horizontal lines show $\chi_{\mathrm{C}}=9\%$ and
$\chi_{\mathrm{H}}=27\%$.

To understand the dependence of the $Q$ rate on the location of the
$\tau$-separators for the four stars, we first fix the number of bins to 4
following Nordlund (1982) and Vögler (2003). Similarly to the experiment
described in section 3.4 of Cernetic et al. (2019), we try all possible
combinations of the $\tau$-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\tau^{\mathrm{ref}}\in[-9.0,0.5]$; for K0V and G2V, 4060 combinations
(grid of 30 equidistant values in $\log\tau^{\mathrm{ref}}\in[-6.5,0.5]$); for
M2V, 4060 combinations (grid of 30 equidistant values in
$\tau^{\mathrm{ref}}\in[-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 $\tau$-separator with maximum optical depth is labelled as
$\tau_{0}$. The indices of the separators increase with height in the
atmosphere ($\tau_{i}>\tau_{i+1}$).

The statistics for each star in our sample show different trends. In general,
combinations of separators that give simultaneously small deviations
$\chi_{\mathrm{C}}$ and $\chi_{\mathrm{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 $\chi_{\mathrm{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 $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ smaller than certain
thresholds indicated individually for each of the stars. For the M2V star the
choice of 4 $\tau$ bins is obviously insufficient. There are only a few
combinations (less than 1%) that lead to $\chi_{\mathrm{C}}<10\%$ and
$\chi_{\mathrm{H}}<20\%$ and there is not a single one that reduces
$\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ below 5% and 20% simultaneously.

How do the deviation measures vary with the location of the three separators?
An illustrative example for the G2V star is shown in Fig. 10, which shows the
dependence of $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ on the location of
the middle separator $\tau_{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\tau_{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
$\tau_{2}$. The position of the highest separator is indicated by the colour,
the black curve corresponding to the smallest value of $\tau_{2}$ and the red
curve to the largest one.

The figure shows that $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ vary
significantly with the location of $\tau_{1}$ if the location of the deepest
separator $\tau_{0}$ is deeper than some critical depth (left-hand columns of
Fig. 10). That critical location of $\tau_{0}$ is located deeper for
$\chi_{\mathrm{C}}$ than for $\chi_{\mathrm{H}}$ ($\log\tau_{0}\simeq-0.5$ for
$\chi_{\mathrm{C}}$ and $\log\tau_{0}\simeq-1.2$ for $\chi_{\mathrm{H}}$). For
$\chi_{\mathrm{C}}$ especially, if $\tau_{0}\simeq-0.5$ the curves flatten
out, meaning that the results are essentially insensitive to the location of
the other two nodes.

If $\tau_{0}$ is deeper than these critical values, the deviation measures are
insensitive to the location of the other two separators only if $\tau_{2}$ is
relatively close to $\tau_{0}$ (the red curve in the left-hand column of Fig.
10). When $\tau_{0}$ is pushed higher up, the deviation measures are
insensitive to the location of the remaining separators and their value
increases (compare two right-hand columns of Fig. 10).

When $\log\tau_{0}\geq-0.5$, the shape of the curves for the
$\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ in Fig. 10 becomes more different.
The cooling deviation shows one minimum around $\log\tau_{1}\approx-2.5$,
while the heating deviation has two minima: one for the locations of
$\tau_{1}$ close to $\tau_{0}$ and another one that moves higher up with
increasing height of $\tau_{0}$. When $\log\tau_{0}=-0.7$, the minimum is at
$\log\tau_{1}\approx-5$, and when $\log\tau_{0}\geq-0.5$, the minimum is at
$\log\tau_{1}\approx-2.5$. Since the location of one of the minima of
$\chi_{\mathrm{H}}$ coincides with the minimum of $\chi_{\mathrm{C}}$, both
can be minimised using the same set of separators even when $\tau_{0}$ is as
deep as $\log\tau_{0}=0.5$.

When the location of the separator $\tau_{0}$ is set above the critical
values, the sensitivity of both $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ to
the location of $\tau_{1}$ and $\tau_{2}$ mostly disappears (see the sequence
from the left to the right panel of Fig. 10). This happens for a higher up
location of $\tau_{0}$ for $\chi_{\mathrm{H}}$ than for $\chi_{\mathrm{C}}$.

Our analysis suggests that, in the case of the G2V star, there are two
strategies for placing the separators so that $\chi_{\mathrm{C}}$ and
$\chi_{\mathrm{H}}$ remain below $4\%$ and $15\%$ respectively. The first
strategy is to set $\log\tau_{0}$ around the continuum formation height,
[-0.5, 0] and $\log\tau_{1}\in[-3,-2.5]$. The second one is to place
$\log\tau_{0}$ so that it includes most of the weak and intermediate
photospheric lines ($\log\tau_{0}\in[-2.4,-1.9]$) and $\log\tau_{1}$ even
higher in $[-5,-3.5]$. In both cases the location of $\log\tau_{2}$ can be
anywhere above $\log\tau_{1}$. The former case (see 2nd and 3rd column)
contains the combinations that produce the minimal $\chi$ 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 $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{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 $\tau_{0}$ located deeper than a certain critical value,
$\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ vary with $\tau_{1}$ similarly to
the G star, showing two minima for $\chi_{\mathrm{H}}$, of which one coincides
with the single minimum of $\chi_{\mathrm{C}}$. There is a little variation of
the deviation curves with $\tau_{2}$ except when it is close to the
$\tau_{0}$. For the F3V star, the optimal combinations that produce
$\chi_{\mathrm{C}}\leq 8\%$ and $\chi_{\mathrm{H}}\leq 15\%$ are (a)
$\log\tau_{0}\in[-0.2,0.5]$ and $\log\tau_{1}\in[-2.5,-2]$ and (b)
$\log\tau_{0}\in[-2.2,-0.7]$ and $\log\tau_{1}\in[-5.5,-4.7]$, with
$\log\tau_{2}$ again anywhere above $\log\tau_{1}$.

For the K0V star, we identified three variants that produce
$\chi_{\mathrm{C}}\leq 5\%$ and $\chi_{\mathrm{H}}\leq 15\%$: (a)
$\log\tau_{0}\in[0,0.3]$ and $\log\tau_{1}\in[-2,-1]$, (b)
$\log\tau_{0}\in[-0.5,-0.2]$ and $\log\tau_{1}\in[-4,-3]$, and
(c)$\log\tau_{0}\in[-1.9,-1.7]$ and $\log\tau_{1}\in[-5,-3]$, with
$\log\tau_{2}<\log\tau_{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\tau_{0}$ there is only one
minimum for the $\chi_{\mathrm{H}}$ varying with $\tau_{1}$ and it does not
coincide anymore with the minimum for $\chi_{\mathrm{C}}$, but with its local
maximum (cf. two left columns in Fig. 11). Secondly, once the $\tau_{0}$
separator is sufficiently high in the atmosphere so that the deviation
measures are insensitive to changes in $\tau_{1}$, $\chi_{\mathrm{C}}$
flattens out as for the other stars but its value quickly increases with the
height of $\tau_{0}$ (see sequence from left to right in the figure), from
$\chi_{\mathrm{C}}\simeq 6\%$ when $\log\tau_{0}=0.2$ to
$\chi_{\mathrm{C}}\simeq 15.5\%$ when $\log\tau_{0}=-0.9$. Higher up, the
value of $\chi_{\mathrm{C}}$ remains larger than 12%, while
$\chi_{\mathrm{H}}$ also starts to increase once it becomes flat (after
$\tau_{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 $\chi_{\mathrm{C}}\leq 9\%$ and
$\chi_{\mathrm{H}}\leq 27\%$ for $\log\tau_{0}\in[0.2,0.5]$ and
$\log\tau_{1}\in[-1.1,-0.6]$, with $\log\tau_{2}<\log\tau_{1}$.

Based on the analysis presented above, we selected optimal four-bins
combinations as reference for other binning strategies that we described in
the following sections: $\log\tau=\left\\{-1.9,-5,-7\right\\}$ for the F3V
star; $\log\tau=\left\\{0,-2.5,-5\right\\}$ for the G2V;
$\log\tau=\left\\{-0.2,-3,-5\right\\}$ for the K0V;
$\log\tau=\left\\{0.3,-0.8,-2.5\right\\}$ for the M2V.

###  5.2 Number of $\tau$-bins

Figure 12: Deviation measures $\chi_{\mathrm{C}}$ (upper panel) and
$\chi_{\mathrm{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, 6, 11, 21, 31, 41, 51, and 101 $\tau$-bins distributed
equidistantly. Circles: deviations for $Q$ computed with the Rosseland
opacity. Pluses: the deviation of the optimal combinations of four $\tau$-bins
(see Sect. 5.1 and Table 4).

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)777In the present and following sections, there
could be combinations of bins where one bin might end up empty. In this case,
when computing the RTE solution, the empty bin is ignored.  that are uniformly
distributed in the relevant range for each of the stars:
$\log\tau^{\mathrm{ref}}\in[-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 Figure 12 $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ are shown as a
function of the number of bins. For F, G, and K – type stars, there is very
little variation in $\chi_{\mathrm{C}}$. The trend is slightly different for
each of the stars individually. While $\chi_{\mathrm{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 $\chi_{\mathrm{C}}$ saturates at around 20 bins. The M star shows
a significant improvement in $\chi_{\mathrm{C}}$ when the number of bins
increases from 4 to 21. For more bins it saturates at $\approx$$5\%$. Although
there is more variation in $\chi_{\mathrm{H}}$ than in $\chi_{\mathrm{C}}$,
the overall trend is similar and, again, the results saturate for more than
$\approx$$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
QODF.

Wray et al. (2006) investigated the performance of the OB and ODF methods in
the context of shock-generated radiation during simulated re-entry of Apollo
AS-501 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}$, $\varkappa$-sorting instead of $\tau$-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 $\chi_{\mathrm{C}}$ and
$\chi_{\mathrm{H}}$ are the minimal deviations that cannot be eliminated by
increasing the number of the $\tau$-separators. Depending on the stellar type,
the minimal $\chi_{\mathrm{C}}$ values span between 2% and 7%, and the minimal
$\chi_{\mathrm{H}}$ span between $\approx$$10\%$ and $\approx$$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
$\chi_{\mathrm{H}}$ value for the F3V star is reduced to only a few per cent
when three separators are set to $\log\tau=\left\\{-1.8,-4.3,-6.6\right\\}$.
This is consistent with locations that produce simultaneous minimum
$\chi_{\mathrm{C}}$ and $\chi_{\mathrm{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. 22). The
$\chi$ differences for this grey solution are shown as coloured circles in
Fig. 12. Their values are larger or similar to the 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.

###  5.3 $\\{\tau,\lambda\\}$-binning

Figure 13: Example of the binning used to compute the $Q$ values for the G2V
star, using different numbers of $\\{\tau,\lambda\\}$-bins. Similarly to Fig.
8, each grey dot corresponds to one ODF substep. For this star, the
$\tau$-separators in the UV in all the experiments are fixed at
$\log\tau=\left\\{-1.8,-4.2\right\\}$. Columns, from left to right, show the
cases for $N$=6,7,8,9 bins. The red numbers show indices of the bins. The case
for $N$=6 has two $\lambda$-separators at 400 and 1600 nm. For more than 6
bins, bins 1 to 5 are kept identical, but more $\lambda$ separators are added
between 400 and 9750 nm. A similar binning set-up is used for the other stars
with modified $\tau$-separators in the UV and the $\lambda$ separator between
the short and long wavelengths for $N$¿6. Figure 14: Deviation measures
$\chi_{\mathrm{C}}$ (upper panel) and $\chi_{\mathrm{H}}$ (lower) versus the
number of bins, for the four stars (colours as in Fig. 1), for the case of
$\\{\tau,\lambda\\}$-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 $\tau$-separators located at the heights of the
bins in the UV (no $\lambda$-separators). Pluses: the deviation of the optimal
combinations of four $\tau$-bins (see Sect. 5.1 and Table 4). Figure 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^{\mathrm{OB}}=\sum_{l}Q_{l}$
for $N=7$. Right panel: comparison of $Q_{8}$ and $Q_{9}$ for $N=9$ against
$Q_{7}$ for $N=7$.  Figure 16: Deviation measures $\chi_{\mathrm{C}}$ (upper
panel) and $\chi_{\mathrm{H}}$ (lower) against the number of bins for the four
stars, computed for the strategy with five fixed $\lambda$-separators at
$\lambda\in(381,726,1383,2636,5023)$ nm and with varying number of
$\tau$-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
$\\{\tau,\lambda\\}$-bins. Circles: deviations for $Q$ computed with the
Rosseland opacity. Pluses: deviations of the optimal combinations of four
$\tau$-bins (see Sect. 5.1 and Table 4).

The OB method using only $\tau$-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 App. 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 $\tau$-binning. Therefore,
$\lambda$-separators between these regions are introduced in addition to the
$\tau$-separators. Such bidimensional distribution of the opacities is non-
uniform in $(\tau,\lambda)$ parameter space: in every wavelength interval the
$\tau$-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 $\\{\tau,\lambda\\}$-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 short-ward and long-ward 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 $\tau$-separators. In addition, the
least opaque bin of the visible + IR opacities is further split with a set
$\lambda$-separators.

We computed the $Q$ values using a $\\{\tau,\lambda\\}$-binning approach
similar to Magic (2014) (see his figure 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 $\lambda$-separator at $\lambda_{0}$. The
value of $\lambda_{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 short-
wavelength opacities are then grouped into three $\tau$-bins, and the long-
wavelength ones into two. The least opaque $\tau$-bin of the long-wavelength
opacities is then split by one or more additional $\lambda$-separators. For
$N=6$ the $\lambda$-separator is placed at 1600 nm. For $N>6$, the
$\lambda$-separators of the least opaque $\tau$-bin are placed equidistantly
in $\log\lambda$ in the range $[\lambda_{0},9570]$ nm. The $\tau$-separators
in the short-wavelength 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 $\tau$-$\lambda$ plane (cf. Figs. 8 and 9).

We computed the $Q$ rate for every combination of $\\{\tau,\lambda\\}$-bins.
The corresponding deviations $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$ are
shown in Fig. 14. The behaviour of $\chi_{\mathrm{C}}$ with the number of bins
is generally similar to that found in Fig. 12. For four stars the
$\chi_{\mathrm{C}}$ deviation decreases with increasing number of bins until
it saturates after around 10–20 bins. The saturation values of
$\chi_{\mathrm{C}}$ are lower than the ones found in Sect. 5.2 (cf. Table 4).
For less than 15 bins ($\Delta\log\lambda>0.1$) $\chi_{\mathrm{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 $\chi_{\mathrm{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 left-hand panel of Fig. 15.
In the right-hand 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 non-linearity 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 $\chi_{\mathrm{C}}$, this leads to a
drop from $\approx$$14\%$ for $N=7$ to $\approx$$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
$\chi_{\mathrm{C}}$.

For the M2V star the $\chi_{\mathrm{H}}$ saturation value is about the same as
in Fig. 12, but for the F, G, and K stars $\chi_{\mathrm{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 $(\tau,\lambda)$ plane. The
number of $\tau$-separators is varied between 1 and 16 and the number of
$\lambda$-separators is fixed at 5 ($\lambda=381,726,1383,2636,5023$ nm). The
$\tau$-separators for all the stars are distributed equidistantly at
$\log\tau^{\mathrm{ref}}\in[-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 $\chi_{\mathrm{H}}$ is larger than the deviation measure
we obtain with the optimised solutions from Sect. 5.1. With 12 bins the values
$\chi_{\mathrm{H}}$ lie between 20 and 30%, similar to those in Fig. 12 with
11 bins distributed only in $\tau$. With finer sampling in $\tau$, the value
of $\chi_{\mathrm{H}}$ is quickly reduced below 10% and saturated at that
level for all four stars.

Table 4: Deviation measures $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{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). In column 5.1 we show the deviations of the 4-bins case with the separators optimally distributed as described in Sect. 5.1. In column 5.2 we show the deviations for the case with 21 bins from Fig. 12. In column 5.3a we show the deviations for the case with 15 bins from Figure 14. Finally, in column 5.3b, we show the deviations for the case with 24 bins from Fig. 16. | 5.1 | 5.2 |  5.3a |  5.3b
---|---|---|---|---
| $\chi_{\mathrm{C}}$ | $\chi_{\mathrm{H}}$ | $\chi_{\mathrm{C}}$ | $\chi_{\mathrm{H}}$ | $\chi_{\mathrm{C}}$ | $\chi_{\mathrm{H}}$ | $\chi_{\mathrm{C}}$ | $\chi_{\mathrm{H}}$
F | 6.8 | 4.8 | 8.0 | 25.4 | 4.7 | 134.1 | 5.4 | 6.0
G | 2.6 | 8.9 | 4.0 | 13.5 | 2.1 | 54.8 | 2.9 | 3.4
K | 2.9 | 5.2 | 2.7 | 6.0 | 1.9 | 39.4 | 1.7 | 2.0
M | 6.0 | 18.9 | 5.9 | 16.1 | 2.2 | 16.7 | 2.2 | 6.3

From this analysis it follows that introduction of the $\lambda$ separators
may reduce the minimum deviation of the OB in particular cases when it is
combined with carefully selected separators in $\tau$. 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 $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{H}}$
computed in the experiments described in this section. The values shown for
the experiments described in Sect. 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 4-bin approach with carefully tuned locations of
the $\tau$-separators (Sect. 5.1) performs similarly or better than the other
approaches for all 4 stars. The only exception is the M star, where taking
into account the wavelength distribution simultaneously with distribution in
$\tau$ for all $\lambda$ 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 Figure 17 $Q$ and $Q/\rho$ 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/\rho$,
since the shape of $Q$ is practically not modified by the relatively slowly
changing density in sub-surface 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/\rho$ show similar information regarding the quality
of the computed $Q^{\mathrm{OB}}$ in respect to $Q^{\mathrm{ODF}}$. Some
differences are expected when choosing $Q/\rho$, 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/\rho$ would be the quantity to
optimise.

Figure 17: Comparison of $Q$ and $Q/\rho$ for the four stars (each row) and
five opacity approaches: ODF (black solid line), grey opacity computed with
Eq. 22 (grey solid line), and the binned opacities for 4 bins (red solid
line), 15 bins (yellow solid line), and 24 bins (black dotted line)
corresponding to columns 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/\rho$ for the same part of the atmosphere. Third column: $Q$ for the higher
layers where the heating occurs. Right column: $Q/\rho$ for the same part of
the atmosphere.

##  6 Conclusions

This work is based on detailed monochromatic opacity tables computed with the
SYNSPEC (Hubeny & Lanz 2011) code using complete up-to-date 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 near-
surface convection as it reduced the number of required radiative transfer
computations by several orders of magnitude while preserving the essence of
the non-grey 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 $\tau$ 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 (App. C). It is possible to find combinations
of small number of bins that significantly reduce the deviations. However,
these combinations are model-dependent and with relatively narrow intervals in
which the bins can be located. Strategies that include preservation of the
wavelength dependence through combined $\tau$\- and $\lambda$-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 $\chi$ values by a factor of 3 (see the bottom row of Table
4), requires an increase in the number of bins, and thus the computing cost by
a factor of 4 or 5. 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 non-solar
metallicities.

An open question left in this work is the importance of the measured
deviations $\chi_{\mathrm{C}}$ and $\chi_{\mathrm{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 section 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 follow-up 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.

###### 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. APG 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. CAP 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. 29. 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.

## References

  * Allende Prieto et al. (2003a) Allende Prieto, C., Hubeny, I., & Lambert, D. L. 2003a, ApJ, 591, 1192
  * Allende Prieto et al. (2022) Allende Prieto, C., Hubeny, I., Lanz, T., & Osorio, Y. 2022, Synple, https://github.com/callendeprieto/synple/blob/master/docs/synple.pdf
  * Allende Prieto et al. (2003b) Allende Prieto, C., Lambert, D. L., Hubeny, I., & Lanz, T. 2003b, ApJS, 147, 363
  * Anders & Grevesse (1989) Anders, E. & Grevesse, N. 1989, Geochim. Cosmochim. Acta., 53, 197
  * Bautista (1997) Bautista, M. A. 1997, A&AS, 122, 167
  * Beeck et al. (2013) Beeck, B., Cameron, R. H., Reiners, A., & Schüssler, M. 2013, A&A, 558, A48
  * Bell (1980) Bell, K. L. 1980, J. Phys. B: Atom. Mol. Phys., 13, 1859
  * Bell & Berrington (1987) Bell, K. L. & Berrington, K. A. 1987, J. Phys. B: Atom. Mol. Phys., 20, 801
  * Borysow et al. (2001) Borysow, A., Jorgensen, U. G., & Fu, Y. 2001, J. Quant. Spec. Radiat. Transf., 68, 235
  * Bruls et al. (1999) Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233
  * Caffau et al. (2008) Caffau, E., Ludwig, H. G., Steffen, M., et al. 2008, A&A, 488, 1031
  * Carbon et al. (1969) Carbon, D., Gingerich, O. J., & Latham, D. W. 1969, in Low-Luminosity Stars, ed. S. S. Kumar, 435
  * Castelli & Kurucz (2003) Castelli, F. & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, Vol. 210, A20
  * Cernetic et al. (2019) Cernetic, M., Shapiro, A. I., Witzke, V., et al. 2019, A&A, 627, A157
  * Collet et al. (2018) Collet, R., Nordlund, Å., Asplund, M., Hayek, W., & Trampedach, R. 2018, MNRAS, 475, 3369
  * Dalgarno (1962) Dalgarno, A. 1962, Geophys. Corp. of America, GCA Tech. Rep., 62
  * Felipe et al. (2010) Felipe, T., Khomenko, E., & Collados, M. 2010, ApJ, 719, 357
  * Fernley et al. (1987) Fernley, J. A., Seaton, M. J., & Taylor, K. T. 1987, J. Phys. B: Atom. Mol. Phys., 20, 6457
  * Gingerich (1969) Gingerich, O., ed. 1969, Theory and observation of normal stellar atmospheres (Cambridge: M.I.T. Press), 377–472
  * Gingerich et al. (1971) Gingerich, O., Noyes, R. W., Kalkofen, W., & Cuny, Y. 1971, Sol. Phys., 18, 347
  * Gustafsson et al. (1975) Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 42, 407
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951
  * Gustafsson & Frommhold (2001) Gustafsson, M. & Frommhold, L. 2001, ApJ, 546, 1168
  * Gustafsson & Frommhold (2003) Gustafsson, M. & Frommhold, L. 2003, A&A, 400, 1161
  * Hayek (2010) Hayek, W. S. 2010, PhD thesis, Australian National University, Canberra
  * Henry (1970) Henry, R. J. W. 1970, ApJ, 161, 1153
  * Hubeny et al. (2021) Hubeny, I., Allende Prieto, C., Osorio, Y., & Lanz, T. 2021, arXiv e-prints, arXiv:2104.02829
  * Hubeny et al. (1994) Hubeny, I., Hummer, D. G., & Lanz, T. 1994, A&A, 282, 151
  * Hubeny & Lanz (2011) Hubeny, I. & Lanz, T. 2011, Synspec: General Spectrum Synthesis Program
  * Hubeny & Lanz (2017a) Hubeny, I. & Lanz, T. 2017a, arXiv e-prints, arXiv:1706.01935
  * Hubeny & Lanz (2017b) Hubeny, I. & Lanz, T. 2017b, arXiv e-prints, arXiv:1706.01937
  * Hubeny & Mihalas (2014) Hubeny, I. & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton University Press)
  * John (1968) John, T. L. 1968, MNRAS, 138, 137
  * Jørgensen et al. (2000) Jørgensen, U. G., Hammer, D., Borysow, A., & Falkesgaard, J. 2000, A&A, 361, 283
  * Khomenko et al. (2017) Khomenko, E., Vitas, N., Collados, M., & de Vicente, A. 2017, A&A, 604, A66
  * Khomenko et al. (2018) Khomenko, E., Vitas, N., Collados, M., & de Vicente, A. 2018, A&A, 618, A87
  * Kurucz (1970) Kurucz, R. L. 1970, SAOSR, 309
  * Kurucz (1993a) Kurucz, R. L. 1993a, Atomic data for opacity calculations., CD-ROM No. 1.
  * Kurucz (1993b) Kurucz, R. L. 1993b, Diatomic Molecular Data for Opacity Calculations., CD-ROM No. 15.
  * Kurucz (1993c) Kurucz, R. L. 1993c, Opacities for Stellar Atmospheres: [+0.0],[+0.5],[+1.0]., CD-ROM No. 2.
  * Kurucz (2002) Kurucz, R. L. 2002, in American Institute of Physics Conference Series, Vol. 636, Atomic and Molecular Data and Their Applications, ed. D. R. Schultz, P. S. Krstic, & F. Ownby, 134–143
  * Kurucz (2014) Kurucz, R. L. 2014, in Determination of Atmospheric Parameters of B-, A-, F- and G-Type Stars, ed. E. Niemczura, B. Smalley, & W. Pych, 39–51
  * Kurucz (2017) Kurucz, R. L. 2017, Can. J. Phys., 95, 825
  * Kurucz et al. (1987) Kurucz, R. L., van Dishoeck, E. F., & Tarafdar, S. P. 1987, ApJ, 322, 992
  * Labs (1951) Labs, D. 1951, ZAp, 29, 199
  * Ludwig (1992) Ludwig, H. G. 1992, PhD thesis, University of Kiel , Germany
  * Ludwig et al. (1994) Ludwig, H. G., Jordan, S., & Steffen, M. 1994, A&A, 284, 105
  * Ludwig & Steffen (2013) Ludwig, H. G. & Steffen, M. 2013, Mem. Soc. Astron. Ital. Suppl., 24, 53
  * Magic (2014) Magic, Z. 2014, PhD thesis, Ludwig-Maximilians University of Munich, Germany
  * Magic et al. (2013) Magic, Z., Collet, R., Asplund, M., et al. 2013, A&A, 557, A26
  * Nahar (1995) Nahar, S. N. 1995, A&A, 293, 967
  * Nordlund (1982) Nordlund, A. 1982, A&A, 107, 1
  * Olson & Kunasz (1987) Olson, G. L. & Kunasz, P. B. 1987, J. Quant. Spec. Radiat. Transf., 38, 325
  * Osorio et al. (2020) Osorio, Y., Allende Prieto, C., Hubeny, I., Mészáros, S., & Shetrone, M. 2020, A&A, 637, A80
  * Palacios et al. (2010) Palacios, A., Gebran, M., Josselin, E., et al. 2010, A&A, 516, A13
  * Peach (1970) Peach, G. 1970, MmRAS, 73, 1
  * Pereira et al. (2013) Pereira, T. M. D., Asplund, M., Collet, R., et al. 2013, A&A, 554, A118
  * Rutten (2003) Rutten, R. J. 2003, Radiative Transfer in Stellar Atmospheres, Utrecht University Lecture Notes
  * Sbordone et al. (2004) Sbordone, L., Bonifacio, P., Castelli, F., & Kurucz, R. L. 2004, Mem. Soc. Astron. Ital. Suppl., 5, 93
  * Skartlien (2000) Skartlien, R. 2000, ApJ, 536, 465
  * Strom & Kurucz (1966) Strom, S. E. & Kurucz, R. L. 1966, J. Quant. Spectrosc. Radiat. Transf., 6, 591
  * Trampedach et al. (2013) Trampedach, R., Asplund, M., Collet, R., Nordlund, Å., & Stein, R. F. 2013, ApJ, 769, 18
  * Vidal et al. (1973) Vidal, C. R., Cooper, J., & Smith, E. W. 1973, ApJS, 25, 37
  * Vögler (2003) Vögler, A. 2003, PhD thesis, University of Göttingen, Germany
  * Vögler et al. (2004) Vögler, A., Bruls, J. H. M. J., & Schüssler, M. 2004, A&A, 421, 741
  * Witzke et al. (2021) Witzke, V., Shapiro, A. I., Cernetic, M., et al. 2021, A&A, 653, A65
  * Wray et al. (2006) Wray, A. A., Ripoll, J. F., & Prabhu, D. 2006, in Studying Turbulence Using Numerical Simulation Database - XI., ed. P. Moin & N. N. Mansour, Proceedings of the 2006 Summer Program, Center for Turbulence Research, Stanford University

##  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
Collision-Induced Absorptions are only taken into account below certain
temperature, in this work chosen to be $8000\mathrm{K}$ (i.e. when molecules
are considered in the EOS).

The free-free cross-sections for H and He+ are computed hydrogenically with
approximated free-free Gaunt factor (see chapter 7.1 of Hubeny & Mihalas
2014). For He, C, N, O, Na, Mg, Al, Si, Ca and Fe the cross-sections are
computed again hydrogenically, but with Gaunt factor of 1. Other free-free
contributions are also taken into account:

  * •

H- from Bell & Berrington (1987).

  * •

H${}_{2}^{+}$ from Gingerich (1969).

  * •

H${}_{2}^{-}$ from Bell (1980).

  * •

He- from the polynomial fit of Carbon et al. (1969) to the data of John
(1968).

The bound-free transitions are mainly computed using hydrogenic cross-section,
data from the Opacity Project TOPbase
888http://cdsweb.u-strasbg.fr/topbase/topbase.html and the Iron Project
999https://cds.unistra.fr/topbase/TheIP.html:

  * •

H-: the cross-section is computed hydrogenically with Gaunt factor of 1 (see
chapter 7.1 of Hubeny & Mihalas 2014).

  * •

H${}_{2}^{+}$: the cross-section is computed using the expression from
Gingerich (1969).

  * •

H: the cross-section for the first nine transitions are computed
hydrogenically with Gaunt factor of 1.

  * •

He: the cross-sections 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 cross-section for 14 transitions is computed hydrogenically with
exact Gaunt factors.

  * •

C: the cross-sections for 104 transitions are interpolated from TOPbase data.

  * •

C+: the cross-sections for 40 transitions are interpolated from TOPbase data.

  * •

N: the cross-sections for 89 transitions are interpolated from TOPbase data.

  * •

N+: the cross-sections for the 51 transitions are interpolated from TOPbase
data.

  * •

O: the cross-sections for 54 transitions are interpolated from TOPbase data.

  * •

O+: the cross-sections for 74 transitions are interpolated from TOPbase data.

  * •

Na: the cross-sections for 32 transitions are interpolated from TOPbase data.

  * •

Na+: the cross-sections for eight transitions are interpolated from TOPbase
data.

  * •

Mg: the cross-sections for 71 transitions are interpolated from TOPbase data.

  * •

Mg+: the cross-sections for the 31 transitions are interpolated from TOPbase
data.

  * •

Al: the cross-sections for 33 transitions are interpolated from TOPbase data.

  * •

Al+: the cross-sections for 81 transitions are interpolated from TOPbase data.

  * •

Si: the cross-sections for 57 transitions are interpolated from TOPbase data.

  * •

Si+: the cross-sections for 46 transitions are interpolated from TOPbase data.

  * •

Ca: the cross-sections for 79 transitions are interpolated from TOPbase data.

  * •

Ca+: the cross-sections for 32 transitions are interpolated from TOPbase data.

  * •

Fe: the cross-sections for 49 transitions are obtained from the Iron Project
data.

  * •

Fe+: the cross-sections for 41 transitions are obtained from the Iron Project
data.

Other contributions to the continuum and bands are also included:

  * •

H, He and H2 Rayleigh scattering from Dalgarno (1962).

  * •

Thomson scattering: $\varkappa=6.65\times
10^{-25}n_{e}/\rho\;\left(\mathrm{cm^{2}g^{-1}}\right)$

  * •

Collision-Induced Absorption of H2-H2 (Borysow et al. 2001), H2-He (Jørgensen
et al. 2000), H2-H (Gustafsson & Frommhold 2003) and H-He (Gustafsson &
Frommhold 2001).

  * •

CH, OH continuous absorption from Kurucz et al. (1987).

Bound-bound 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 line-to-continuum opacity ratio). These are mainly from
Kurucz line lists 101010http://kurucz.harvard.edu/linelists.html, updated by
data from the National Institute of Standards and Technology (NIST) when
available. SYNSPEC also uses data from EXOMOL 111111https://www.exomol.com for
some of the molecular transitions. In the case of the molecules, around a
total of $2\times 10^{7}$ transitions for H2, CH, NH, OH, NaH, MgH, SiH, CaH,
CrH, FeH, C2, CN, CO, MgO, AlO, SiO, CaO, VO, H2O and TiO are included (see
Fig. 5). In the case of the atoms, SYNSPEC includes a line list with around
$2\times 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 2 ions of Be, Nb, Rh, Zr, and Ru.

  * •

Neutral and first 3 ions of B, C, and Mo.

  * •

Neutral and first 5 ions of N, O, F, Ne, Na, Mg, Al, Si, P, S, Cl, Ar, and K.

  * •

Neutral and first 7 ions of Zn.

  * •

Neutral and first 8 ions of Ca, Sc, Ti, V, Cr, Mn, and Co.

  * •

Neutral and first 9 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).

Figure 18: 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
$\lambda\in[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. Figure 19: Same as Fig. 18, for
the step defined by $\lambda\in[650,660]$ nm containing hydrogen H$\alpha$.
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 3 and 4) for the
three hotter stellar types. Figure 20: Same as Fig. 18 for the step defined by
$\lambda\in[580,590]$ nm containing NaI D1 and D2 line. Figure 21: Same as
Fig. 18 for the step defined by $\lambda\in[193,197.5]$ nm. The differences in
the UV between the two tables are large in all the substeps and for all
stellar models. Figure 22: Radiative energy exchange rate $Q_{i}$ computed
from the four model atmospheres using the opacity shown in Fig. 18 and
integrated over the interval $\lambda\in[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). Figure 23: Same as Fig. 22, but for the
$\lambda\in[650,660]$ nm interval that contains the hydrogen H$\alpha$ line.
Computed using the opacity shown in Fig. 19. Figure 24: Same as Fig. 22, but
for $\lambda\in[580,590]$ nm, interval that contains the NaI D1 and D2 line.
Computed using the opacity shown in Fig. 20. Figure 25: Same as Fig. 22, but
for $\lambda\in[193,197.5]$ nm, computed using the opacity shown in Fig. 21.

##  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 $\tau$-bins (Vögler et al. 2004). To explain why this is the case,
we use the checkerboard distribution of the bins introduced in Fig. 26,
similar to that used in Sect. 5.3, but now with 4 bins in $\tau$ and 6 in
$\lambda$. The number of non-empty $\tau$-bins changes depending on the
$\lambda$-bin. For example, in the case of the F3V star there are 4
$\tau$-bins for $\lambda<300$ nm, while there are 2 $\tau$-bins for
$\lambda\in[1000,2000]$ nm (Fig. 26). For each star, in Fig. 27 we compare the
binned opacity of pairs of wavelength ranges selected to better illustrate
this behaviour. We compare $\lambda$-ranges that have the same number and
location of non-empty $\tau$-bins. For each star and any of the $\tau$-bins,
the opacity from one of the two $\lambda$-ranges vary with height in a
different way than the opacity from the other. Even for arbitrarily increased
number of bins, binning only in $\tau$ 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 $\lambda$, since each contributor to the opacity may
have its own variation with height.

Figure 26: 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
($\lambda<300$, $\in[300,500]$, $\in[500,700]$, $\in[700,1000]$,
$\in[1000-2000]$, and $>2000\,\mathrm{nm}$). All values within one wavelength
range are shown in the same colour. The $\tau$-separators are selected from
the optimal combinations found in Sect. 5.1 for each star:
$\log\tau=\left\\{-1.9,-5,-7\right\\}$ for the F3V star;
$\log\tau=\left\\{0,-2.5,-5\right\\}$ for the G2V;
$\log\tau=\left\\{-0.2,-3,-5\right\\}$ for the K0V;
$\log\tau=\left\\{0.3,-0.8,-2.5\right\\}$ for the M2V. Figure 27: Opacity
$\overline{\varkappa}_{l}$ of each $\tau$-bin in pairs of $\lambda$-ranges for
the four stars (spectral type at the top left of each panel) versus the
geometrical hegiht (bottom axis of each panel) and optical depth (top axis).
The $\lambda$-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 $\lambda$-ranges in Fig. 26, with the values of each
range in the legend.

Figure 28 shows how the opacity computed with the OB method converges to the
Rosseland opacity in each $\tau$-bin,
$\overline{\varkappa}_{l}=\overline{\varkappa}^{\mathrm{Ro}}_{l}$ (dotted
black lines), in the limit of large enough optical depths. We show only the
case for the F3V star and the range for $\lambda<300\,\mathrm{nm}$ to avoid
overcrowding the plot. However, the conclusions are valid for any
$\lambda$-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 28 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}_{\mathrm{Ro}}=\int\frac{\partial\mathrm{B}}{\partial\mathrm{T}}d\lambda\Bigg{/}\left(\int\frac{\partial\mathrm{B}}{\partial\mathrm{T}}\frac{1}{\varkappa_{\lambda}}d\lambda\right),$
(22)

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 29 shows that the combination of the opacity bins in the
Rosseland fashion,

$\left.\sum_{l}\,\left.\frac{\partial B}{\partial
T}\right|_{l}\right/\left(\sum_{l}\,\left.\frac{\partial B}{\partial
T}\right|_{l}\,\frac{1}{\overline{\varkappa_{l}}}\right),$ (23)

yields the correct value for the Rosseland mean for deep enough optical
depths. For these optical depths (see Fig. 28 and Eq. 20)

$\overline{\varkappa}_{l}\approx\overline{\varkappa}^{\mathrm{Ro}}_{l}=\left.\frac{\partial\mathrm{B}}{\partial\mathrm{T}}\right|_{l}\Bigg{/}\left(\sum_{i(l)}\Delta\lambda_{i}\frac{\partial\mathrm{B}_{i}}{\partial\mathrm{T}}\sum_{j(i,l)}\frac{\omega_{j}}{\varkappa_{i,j}}\right),$
(24)

thus, the combination of the bins following Eq. 23 gives

$\left.\sum_{l}\,\left.\frac{\partial B}{\partial
T}\right|_{l}\right/\left(\sum_{l}\,\left.\frac{\partial B}{\partial
T}\right|_{l}\,\frac{1}{\overline{\varkappa_{l}}}\right)=\\\
\left.\sum_{l}\,\sum_{i(l)}\Delta\lambda_{i}\frac{\partial\mathrm{B}_{i}}{\partial\mathrm{T}}\right/\left(\sum_{l}\sum_{i(l)}\frac{\partial\mathrm{B}_{i}}{\partial\mathrm{T}}\sum_{j(i,l)}\frac{\omega_{j}}{\varkappa_{i,j}}\right).$
(25)

This equation is similar to Eq. 22, 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

$\left.\sum_{l}\,\left.\frac{\partial B}{\partial
T}\right|_{l}\right/\left(\sum_{l}\,\left.\frac{\partial B}{\partial
T}\right|_{l}\,\frac{1}{\overline{\varkappa_{l}}}\right)\equiv\\\
\left.\sum_{l}\,\int_{\lambda(l)}\frac{\partial\mathrm{B}_{\lambda}}{\partial\mathrm{T}}d\lambda\right/\left(\sum_{l}\int_{\lambda(l)}\frac{\partial\mathrm{B}_{\lambda}}{\partial\mathrm{T}}\frac{1}{\varkappa_{\lambda(l)}}d\lambda\right),$
(26)

where $\lambda(l)$ corresponds to the wavelengths that enter the $l$th bin.
Binning in $\tau$, $\lambda$, or both always is equivalent to group only
separating in $\lambda$, since each opacity point corresponds to a different
wavelength121212In the case of the ODF, the substeps lose their wavelength
dependence, but the ODF formalism takes care of this.. This property as well
as understanding the integrals of the opacity in terms of Riemann or Lebesque
integration explains the result shown in Fig. 29. Similarly, the Planck mean
should be recovered for higher layers.

Figure 28: Opacity $\overline{\varkappa}_{l}$ for each $\tau$-bin in
$\lambda<300\,\mathrm{nm}$ computed with Eq. 20 (solid blue curves) compared
to the opacity
$\overline{\varkappa}_{l}=\overline{\varkappa}^{\mathrm{Ro}}_{l}$ (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. 22. The bottom and top axis
show the geometrical height and optical depth, respectively.  Figure 29:
Rosseland opacity mean computed with Eq. 22 (black dashed line) and Eq. 23
(red thin solid line) for the F3V star. The bottom and top axis show the
geometrical height and optical depth, respectively.