© The Authors 2023

1 Introduction

Radiative transfer is one of the key factors governing the structure and dynamics of stellar atmospheres and it has to be taken into account with special care when constructing a model atmosphere (Hubeny & Mihalas 2014). Detailed interactions between radiation and matter are quantitatively described by the opacity, that is by the absorption coefficients of all contributing (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 ≈10⁹ lines and ≈10⁶ 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 thermody-namic variables (for example, temperature T and mass density ρ); it can therefore be pre-computed 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 pre-computed without the contribution of the non-LTE (NLTE) lines, which are computed explicitly. The opacities are computed as byproducts in the codes and projects for modelling stellar atmospheres: the ATLAS family of codes (Kurucz 1970; Castelli & Kurucz 2003; Sbordone et al. 2004; Witzke et al. 2021), MARCS (Gustafsson et al. 1975, 2008), PHOENIX (Gustafsson et al. 1975, 2008), SYNSPEC/TLUSTY (Hubeny & Lanz 2011; Hubeny et al. 2021), and others.

Second, in the energy conservation equation of the system of magnetohydrodynamic equations, the radiative energy exchange is accounted for by a wavelength-integrated source term Q. The distribution of radiative heating and cooling with wavelength is irrelevant, at least as long as the one-fluid magnetohydrody-namics (MHD) approach is used. Therefore, instead of solving monochromatic RTE and then integrating the results, an approximate solution may be obtained by first combining or integrating the opacities in wavelength ranges and then solving the RTE for a number of representative opacity values. This idea is implemented in the opacity distribution function (ODF) defined by Strom & Kurucz (1966). In their approach the monochromatic opacities are divided into hundreds of wavelength segments. The opacities in each segment are then sorted by magnitude to obtain a continuous monotonie distribution that is then discretised by averaging over several (around 10) substeps. In that way the RTE can be solved for several thousands of ODF values instead of for millions of monochromatic opacities.

An important underlying assumption in the ODF approach is that there are no substantial wavelength shifts in the spectrum along the different heights in the atmosphere. There are situations in which this assumption fails (see chapter 17.6 in Hubeny & Mihalas 2014), especially for plasmas with velocity fields. Nevertheless, even in moving media, if we are not interested in detailed spectra and if the velocities are not larger than the thermal velocity, we can still solve the RTE using ODFs.

Pre-computed 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 speed-up 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 source¹ 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 prep.) with the MANCHA code (Khomenko et al. 2017, 2018), flexibility in selecting opacity contributors is of particular importance for us.

In Sect. 2 we briefly introduce a set of the ID models that represent our stars and specify the radiative transfer details. Computation of the monochromatic opacities using the SYNSPEC code are described in Sect. 3. In Sect. 4 the essentials of the ODF method are summarised: ODFs are constructed from the monochromatic opacities computed with SYNSPEC, and the results are compared with the results of the ATLAS code. The same data set is used to study strategies for OB in Sect. 5: we first test the importance of the location of the bin separators in optical depth (Sect. 5.1), then we explore the influence of the number of bins in optical depth (Sect. 5.2), and, finally, we study two different strategies for combined binning in optical depth and wavelength (Sect. 5.3). Our conclusions are presented in Sect. 6.

Fig. 1 Temperature stratification of the F3V (green), G2V (yellow). K0V (orange), and M2V (dark red) 1D hydrostatic models used in this work. The same colour coding is used throughout the paper to label results for the four stellar types. The grid used to solve the RTE is marked by dots.

2 Method

2.1 Model atmospheres

To test the OB for a series of main sequence spectral types, a set of 1D stellar atmospheric models is used. These atmospheres are computed assuming hydrostatic equilibrium and using the mixing length theory to account for the convective energy transfer. The models are initiated with radiation in the grey approximation and corrected to have representative temperature gradient at the top of the atmosphere by scaling the Hopf function q(τ) (see chapter 17 from Hubeny & Mihalas 2014, and references 507– 510 there) of the Harvard-Smithsonian Reference Atmosphere (Gingerich et al. 1971) with the effective temperature of each star, following the expression . 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⁻²) 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

(1)

where I_v(z, ν) is the specific intensity (erg s⁻¹ cm⁻² Hz⁻¹ ster⁻¹ ) at frequency ν in the direction µ and at the geometrical height z in the atmosphere; B_v is the Planck function; ρ(z) is the mass density; and ϰ_ν (z) is the monochromatic absorption coefficient per unit mass (cm² g⁻¹ ; 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),

(2)

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

We solve the RTE in the vertical direction (µ = ±1). At each frequency ν, the formal solution for outward and inward intensity (I⁺, I⁻) in LTE is:

(3)

where i is the height index for each point in the atmosphere, being zero at the bottom and increasing upwards; is the optical depth at the height with index i (dτ_v(z) = ρ(z)ϰ_ν (z)dz); and and . To solve the RTE, we apply the short characteristics method (Olson & Kunasz 1987), using the linear approximation of the Planck function:

(4)

For each frequency ν we compute the coefficients for the local and the upwind points

(5)

Since the method is formulated in a way that it is symmetrical along a ray direction, the coefficients are the same for the upward and inward intensities, and the only change needed is the specific Δτ_ν in each case ( for the upward intensities; 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_v (erg s⁻¹ cm⁻² Hz⁻¹ ster⁻¹) and the flux ℱ_v (erg s⁻¹ cm⁻² Hz⁻¹) can then be calculated for a vertical ray from the known I_v:

(6)

(7)

Finally, the monochromatic radiative energy exchange rate Q_v (erg cm⁻³ s⁻¹ Hz⁻¹), which accounts for the radiation sources and sinks in the energy conservation equation, can be computed either as the divergence of the radiative flux,

(8)

or directly from the mean intensity,

(9)

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

(10)

where τ_h = 0.1. As it is explained in Vögler (2003), this smooth transition between and avoids the accuracy problem that has in the optically thick regime (where J_v approaches B_v) and the errors coming from small variations of the orientation of the flux in the optically thin part of the atmosphere (where flux should be nearly constant), that are enhanced by the gradient in .

3 Monochromatic opacity

For the computation of the monochromatic opacity ϰ_v, we use synple , a Python wrapper for SYNSPEC . SYNSPEC is a general spectrum synthesis code that has been used to solve the RTE in different astrophysical scenarios, including cool stars (see, for example, Allende Prieto et al. 2003b,a and Palacios et al. 2010; Osorio et al. 2020). The main purpose of SYNSPEC is to synthesise spectra for a given model atmosphere, but it can also be used to compute lookup tables of monochromatic opacities for a given grid of thermodynamic quantities. The code solves the equation of state (EOS), and uses the number densities to compute the continuum and line opacities. The EOS from SYNSPEC is fully consistent with the results of our own EOS solver (Vitas et al., in prep.), which is used to produce the EOS tables for the 3D simulations of stellar atmospheres with the MANCHA3D code (Felipe et al. 2010; Khomenko et al. 2017, 2018; Modestov et al., in prep.).

The EOS in SYNSPEC is solved for a specified grid of temperature T and mass density ρ. The number of nuclei per species (fixed by the abundances) and charge are conserved. The code solves the chemical equilibrium equations to determine the populations for 38 atomic species. For temperatures lower than a certain value (in this work chosen to be 8000 K), neutral atoms and singly ionised atoms are included in the EOS (not taking into account higher ionisations) and the code also considers molecular formation (including 503 molecular species). For higher temperatures, molecules are not included in the EOS, and higher ionisations are computed for the atoms. For more details on the solution of the EOS in SYNSPEC (and other details and the general use of the code), we refer the reader to Section 2.6 from Hubeny et al. (2021), and Hubeny & Lanz (2Ol7a,b).

Once the EOS is solved, the opacity is computed in the same (T,ρ) grid, using the atomic and molecular data from several sources (see Appendix A). The data used by the code can be customised by the user by changing, for example, the line lists for the 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 Appendix A, we computed a lookup table using synple, with a microturbulence velocity of 2 km s⁻¹ and for the solar abundances from Anders & Grevesse (1989). These solar abundances are outdated but preferred in our study since they were used in the computation of the reference ODFs to which they are subsequently compared. The T-axis of the grid has 60 values in the range Τ ϵ [1995, 125 900] K, with constant step Δ log T ≃ 0.03 Κ (in the present text, log ≡ log ₁₀). The ρ-axis has 35 values in the range ρ ∈ [10⁻¹³, 3.2 × 10⁻³] g cm⁻³ with Δ log ρ ≃ 0.3. The table is computed for wavelengths in the range λ ∈ [20, 9.5 × 10⁴] nm, with Δ log λ = 10⁻⁶ (λ in nm). In total there is around 7.7 × 10⁹ data points in the table, stored in an HDF5 file of 230 GB. Another table is computed for the same wavelength grid but with finer temperature and density axes (same number of points, but reduced ranges, Τ ∈ [1995, 19 900] Κ and ρ ∈ [10⁻⁹,3.2 × 10⁻³] g cm⁻³). This table is more suitable for the M2V star, thus reducing the interpolation errors in the RTE solver.

The total opacity (continuum and lines) in the wavelength range from 100 to 2500 nm is shown in Fig. 2. It is interpolated for the temperatures and mass densities of the four 1D stellar models. The opacity clearly changes with different heights in the stellar models and with spectral type. To check where the radiation at each wavelength contributes to the emergent radiation, we mark the geometrical heights at which τ_ν = 1 (thick line). This curve includes the contribution from both the line and the continuum opacity. Since the wavelength range is wide, we divided it into 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 . The ionisation edges and the typical shape of the bound-bound and bound-free H⁻ opacity coefficient are visible for the curve. In the UV the line forest elevates the height of the τ_λ = 1 curve. The strong spectral lines, especially the hydrogen series, are clearly visible as well as their relative decrease in strength with decreasing effective temperature. In the M star, the τ_λ = 1 curve is shaped by millions of molecular transitions.

Following the method from Sect. 2, we solve the RTE to obtain Q_v. Once Q_v is computed, we can integrate it in frequency using the trapezoidal rule:

(11)

The computed Q_v is shown in Fig. 3, for the wavelength range and set of stellar models as in Fig. 2. The radiative energy exchange happens mostly around the stellar surface. For decreasing effective temperature, the absolute value of Q_v is large for a wider wavelength range and the significance of the contribution to Q_v shifts from the visible for the F star towards the IR for the M star. Figure 3 also shows how the relative magnitude and area covered by significant Q_v for the heating in the atmosphere (Q_v > 0) and the cooling (Q_v < 0) changes for the different spectral types. For the F star the maximum heating is around an order of magnitude smaller than the magnitude of the most intense cooling, and peaks in the 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 pre-computed 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⁻³ s⁻¹) and not the ratio Q/ρ (erg g⁻¹ s⁻¹), which would put more emphasis on the radiation exchange at lower densities.

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

Fig. 3 Monochromatic radiative energy exchange rate Qv (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 |Qv| to max |Qv|. Blue tones in the colour-map correspond to negative Qv, i.e. cooling in the atmosphere, and red tones correspond to positive Qv, i.e. heating.

4 Opacity distribution functions

4.1 Construction of the opacity distribution function

The ODF method (Labs 1951; Strom & Kurucz 1966) was invented to allow the accurate calculation of the radiative losses due to the millions of opacity contributors in a computationally efficient way. In the 1960s the solution for the RTE was prohibitive for more than a few wavelengths in 1D model atmospheres. Despite the exponential growth in computing power since then, the same problem of efficiently and accurately accounting for the millions spectral lines is still present today when the RTE is solved in 3D 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 λ ∈ [λ₂, λ₃], is shown at the left panel of Fig. 4). Then, all the opacities within each of the steps are sorted by increasing magnitude (see the middle panel of Fig. 4). At this point the direct mapping between the monochromatic opacity values and their wavelengths is lost. To discretise this monotonic distribution within each step, a number of substeps are defined specifying a set of dimensionless weights ω_j(Σ_j ω_j = 1). The length of a substep Δλ_i,j is calculated in terms of the length of the step Δλ_i and the weight ω_j, so that Δλ_i,j = ω_jΔλ_i ². Finally, the opacities are averaged computing the arithmetic mean within each substep (see the right panel of Fig. 4).

This procedure can be performed for any set of monochromatic opacities computed for a pair of thermodynamic quantities. One of them is usually the temperature T and the other is either the density ρ (as in SYNSPEC³) or the gas pressure P (as in the ATLAS code, Kurucz 1970). Kurucz (1993c) created ODF tables⁴ in the described way for specific grids of (T, P). We apply the same strategy to the monochromatic opacities computed using SYNSPEC in a (T, ρ) grid, for the two monochromatic opacity tables described (see Sect. 3). The wavelength range of our SYNSPEC computation and derived ODF is shorter than the range covered by the ATLAS ODF, λ ∈ [9, 10⁷] nm. When the ATLAS table is used to compute the bolometric RT quantities in the different model atmospheres, there is no difference between using it over its entire λ range and in the reduced range of our SYNSPEC computations. This is not surprising since the emission of cool stars peaks in spectral regions far from the extremes that we left out from the full wavelength interval of the ATLAS tables. For our ODF calculation, the wavelength range is divided into 291 steps, which have exactly the same locations and the same 12 weights for the substeps as ATLAS. The length of the steps is approximately proportional to the wavelength (see the table from Kurucz 1993c, to check the exact values); and the substeps are constructed using the weights ω_j = {0.1, 0.1,0.1, 0.1,0.1,0.1,0.1,0.1, 0.1,0.05,2/60,1 /60}. This kind of 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 cm² g⁻¹, just as it is the case for ATLAS.

The tables from Kurucz (1993c) have been widely used in the context of MHD simulations of solar and stellar atmospheres (see, for example, Vögler 2003; Khomenko et al. 2018; Beeck et al. 2013, and many others). The atomic and molecular line lists used to compute these tables (Kurucz 1993a,b) are sufficiently complete and accurate to make these simulations closely resembling various observational diagnostics. However, apart from the molecules that are included (H₂, HD, MgH, NH, CH, SiH, OH, CH, C₂, CN, CO, SiO), many are missing (e.g. VO and TiO) which are dominant opacity sources in the spectra of M stars. SYNSPEC, however, allows for a more complete set of molecules (see Appendix A and Fig. 5) and includes an 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.

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

4.2 Comparison of opacity distribution function from SYNSPEC and ATLAS

While there is overall good match between the two data sets, there are also significant differences. To understand better these differences, we compare the opacity of each step and substep of the ODFs interpolated for the values of T and ρ along the four 1D model of atmospheres. Four examples representing typical differences are shown in Figs. B.1–B.4 (where black lines represent the opacity from ATLAS along the models, and the coloured lines on top of them represent the opacity from SYNSPEC).

In parallel, we also compare the bolometric Q rate computed using both ODFs. We apply the method presented in Sect. 2 to solve the RTE and compute Q_i,j for every step i and substep j; we then integrate in frequency following the ODF formalism:

(12)

where Q_i is the rate integrated in wavelength in the ith step.

Comparing the values of Q computed from different data sets is not straightforward. An example of Q computed for the model of the M star using two opacity sets and their absolute difference is shown in Fig. 6. In the figure we identify the typical negative feature that corresponds to the strong cooling at the bottom of the photosphere and the relatively smaller positive one that shows the radiative heating (blanketing effect). Higher up and deeper down in the atmosphere, the values of Q can be several orders of magnitude smaller. While the radiative losses at these heights may still be significant, we focus here on the dominant components around the optical depth unity. Wherever the values of Q are small, the relative difference computed between them may be exaggerated with respect to the importance of these differences from the energy balance in the atmosphere. Moreover, as the sign of Q flips in the region of interest (the interval of heights around τ = 1) and the two Q values computed from ATLAS and SYNSPEC do not necessary flip the sign at exactly the same height, the relative difference around that height might also be misleadingly high. Therefore, instead of computing the relative differences in the entire domain, we focus on the two dominant features and measure the deviation for the cooling and heating part separately.

For each of them we compare the area A of the difference of Q computed for each of the ODF data sets normalised to the area of one of them:

(13)

(14)

where we choose as Q⁽¹⁾ and Q⁽²⁾ 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 and . The height at z_C→H is the height where Q changes sign for the first time above the cooling component; z_b and z_t are, respectively, the highest point in the atmosphere under the surface (z < 0 in the plots) and the lowest point over the surface (z > 0) where |Q| < 2 × 10⁻⁴ |min(Q)| (see Fig. 6). If A (Q)_H < 1% A (Q)_C, χ_H is not computed (these cases are marked as χ_H = - - - % in the figures).

The radiative energy exchange terms Q_i computed from the four models using the same steps as in Figs. B.1–B.4 are shown in Figs. B.5–B.8. The corresponding values of the deviation measures are indicated in each panel.

When comparing the opacity from the ODF of the two data sets, we see a general match in behaviour and values in most of the visible and infrared intervals for the F, G, and K stars. An example for a step with a good match (λ ∈ [590, 600] nm) is shown in Fig. B.1. The Q_i values computed in the same step closely overlap for these stars with all χ_C and χ_Η between 1.4% and 3.1% (see, for example, Fig. B.5). For the 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α (Fig. B.2) or NaI D1 and D2 (Fig. B.3). Similar discrepancies are found in the steps containing other strong lines like the CaII H and K lines, the G band, HI β, MgI b1, the CaII triplet, and the most intense lines from the Paschen, Brackett, and Pfund series. In the case of all the hydrogen series the opacity from ATLAS is higher than the one from SYNSPEC in the deeper part of the atmospheres, but this reverses in the upper part (see Fig. B.2). For these hydrogen lines, such opacity behaviour is present in all the substeps, although the reversal occurs deeper in the less opaque substeps. We notice a similar tendency in the G band. Another trend can be seen for the steps containing CaII lines, MgI b1 and NaI D1 and D2, for which the opacity from SYNSPEC is higher than the one from ATLAS at all the heights of the atmospheres for the most opaque substeps (see Fig. B.3). The number of affected substeps varies between the two and four most opaque, depending on which intense line falls within the step. The mismatch is affected by the number of strong lines contained in the interval, the magnitude of their opacity with respect to the continuum, and their width relative to the length of the interval. In the less opaque substeps, the two data sets give very close results.

In the steps that contain strong spectral features, the radiative energy exchange rates generally show χ_C,χ_H ≲ 5–10% for the F, G, and K stars (see, for example, Figs. B.6, B.7). An exception happens for the G star for the step containing hydrogen Hα where χ_Η ≃ 25% (see Fig. B.6). However, in the latter example, the area of the heating feature is only 0.02 times the area of the cooling, making this discrepancy less important.

The differences between SYNSPEC and ATLAS for the strong atomic lines could be, for example, due to the wavelength resolution in the monochromatic opacity tables or the broadening parameters used in the line synthesis (e.g. Van der Waals broadening or damping constants). In the case of the hydrogen lines, the difference is explained by the treatment of the Stark broadening. For this project, we used approximate Stark broadening profiles after Hubeny et al. (1994), while in ATLAS the tables from Vidal et al. (1973) are used. Without having access to the monochromatic opacities from ATLAS it is difficult to pinpoint the exact reason of the discrepancies for the rest of the species. Although both codes are open source, their complexity makes it very difficult to find the exact source of differences by comparing the actual source code.

Apart from the spectral features discussed above, the comparison of the two opacity sets in the visible and IR for the F, G, and K stars reveals an excellent agreement. At the same time, there is a notable mismatch for the opacity in the UV between the two data sets (see an example in Fig. B.4). These discrepancies are also evident in the corresponding values of Q_i (see an example in Fig. B.8), for which we find differences of around χ_C,χ_H ≃ 5–40% for λ ≳ 190 nm (as in the example from Fig. B.8), and around χ_C,χ_H ≃ 30–80% for λ < 190 nm (in other steps in the ultraviolet).

The differences between the two data sets in the UV are likely due to the use of photoionisation cross-sections from the Opacity Project⁵ 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 χ_C and χ_Η are also larger, for most of the steps in the visible and NIR they are in the range ≃10–35% (≃40–60% in the worst cases).

As previously discussed, the discrepancies between the two tables for the 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 χ_Η ≃ 77%, in contrast to the values χ_C < 3% and χ_Η < 10% found for the F, G, and K stars. To demonstrate the effect of certain molecules on the discrepancies between the data sets, we recompute the SYNSPEC table and Q, excluding TiO and VO from the line lists for the M star (blue curve in the bottom right panel of Fig. 7). In this case, the value of χ_H is considerably reduced, but it is still relatively high (≃29%). Even if the selection of molecules in the two codes is identical, differences between the opacities are expected owing to different line lists and other atomic and molecular data. For example, the use in SYNSPEC of more up to date data for the opacities, in particular for H₂O and TiO, from EXOMOL⁶, or different partition functions and oscillator strengths.

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

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

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

5 Opacity binning

A further approximation of the opacity computation is the OB method initially introduced by Nordlund (1982). Its implementation in the MURaM and CO⁵ BOLD codes is described by Vögler (2003) and Ludwig (1992) respectively. In this method the opacities are sorted taking into account their distribution with the discretised optical depth τ^ref of a representative 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 τ-separators (). In a variation of the method (see Magic 2014) the wavelength dependence is partially preserved by grouping opacities by both optical depth and wavelength. In that case each bin is defined by two pairs of separators, one in optical depth (τ-separators) and one in wavelength (λ-separators).

We start from the ODF constructed as it is described in Sect. 4. The optical depth is computed for each step i and sub-step j of the ODF for all geometrical heights z in a model atmosphere as

(15)

where ρ refers to the mass density in the model atmosphere (and not the density in the grid of the opacity table) and the height dependence in ϰ_i,j(z) denotes that the opacity from the ODF is interpolated into the temperatures and densities from the model atmosphere. The opacities are then grouped so that opacity of the ODF step i and substep j belong to bin l if , where z_i,j is the height at which τ_i,j = 1. Additionally, when the λ-separators are used, every ODF value is sorted by wavelength so that the central wavelength of an ODF step falls between a pair of λ-separators. Figure 8 shows examples of such grouping of the ODF values into five bins for our four stellar models. There are two λ-separators approximately dividing the wavelength axis into ultraviolet, visible, and infrared and there is one τ-separator for the UV (log τ = −3) and a joint one for the visible and IR (log τ = −1). The corresponding 2D histograms of the ODF points are shown in Fig. 9. The minimal optical depth in the reference models vary with the stellar type and wavelength. We assign the top geometrical height of the atmosphere z_i,j = z_top to all ODF values with τ_i,j(z_top) > 1.

After the opacities of the ODF have been distributed into bins, for every bin l (that includes the ODF steps i = i(l) and substeps j = j(i,l) contained in the bin) and for every T, ρ of the grid of the ODF we compute the Planck function,

(16)

its derivative,

(17)

the Planck mean opacity,

(18)

and the Rosseland mean opacity,

(19)

where Δλ_i, ω_j, and ϰ_i,j are, respectively, the length of the wavelength step, the weight of the substep, and the opacity of the ith step and jth substep of the ODF. B_i is the Plank function in the middle of the wavelength step for each temperature T in the (T, ρ) grid used to compute the ODF.

Finally, following Ludwig (1992), and are combined to get a mean value of the opacity for each bin:

(20)

where the transition between thin and thick optical regimes is at the threshold optical depth τ_thr = 0.35 and τ_l is computed using pressure p in the ODF grid and a stellar surface gravity g, approximately as . With this definition, the opacity from each bin converges to the grey approximation for large enough optical depths (see the discussion about Fig. C.3).

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

(21)

where Q_l refers to partial rate of the bin l computed using Eq. (10).

We apply the described procedure to the ODF computed with SYNSPEC monochromatic opacities presented in Sect. 4. There are several important free parameters in the procedure: the model atmosphere itself, the number of bins in τ and λ, and the locations of separators between the bins. There is no obviously intuitive choice for these parameters. We designed numerical experiments to test various alternatives. As a deviation measure of the radiative energy exchange rate Q⁽²⁾ (computed using the binned opacity) with respect to the reference rate Q⁽¹⁾ (computed with ODF), Vögler et al. (2004) plotted the height-dependent absolute difference of Q⁽¹⁾ and Q⁽²⁾ divided by the density stratification ρ of the atmosphere. The scaling with density is important to compare the radiative rates with other energy sources and sinks in the simulations. Alternatively, Magic (2014) opted for computing max |Q⁽¹⁾ − Q⁽²⁾|/ max |Q⁽¹⁾|, 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⁽¹⁾) and the values that are tested (Q⁽²⁾) 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 multiparameter study we break the problem down by first studying the influence of the location of the τ-separators for a fixed number of bins (4) and with no bins in the wavelength direction. We then allow the number of τ bins to change. Finally, we test several cases including the λ bins.

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

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

5.1 τ-bin location

To understand the dependence of the Q rate on the location of the τ-separators for the four stars, we first fix the number of bins to four following Nordlund (1982) and Vögler (2003). Similarly to the experiment described in Sect. 3.4 of Cernetic et al. (2019), we try all possible combinations of the τ-separators from a discretised grid of optical depths for each of the stars. For the F3V star this gives us 9880 combinations selected from a grid of 40 equidistant values in log τ^ref ∈ [−9.0, 0.5]; for K0V and G2V, 4060 combinations (grid of 30 equidistant values in log τ^ref ∈ [−6.5, 0.5]); for M2V, 4060 combinations (grid of 30 equidistant values in τ^ref ∈ [−4.5, 0.5]). The grid is customised for each star to cover properly the relevant range in optical depths (see the vertical axis of Fig. 8). The τ-separator with maximum optical depth is labelled as τ₀. The indices of the separators increase with height in the atmosphere (τ_i > τ_i+1).

The statistics for each star in our sample show different trends. In general, combinations of separators that give simultaneously small deviations χ^C and χ^H are easy to find for the G and K stars, while they become less frequent in the case of the F star and quite rare for the M star. The number of cases leading to the deviations below certain threshold is given in Tables 1 and 2 separately for the cooling and the heating, and in Table 3 when both components are constrained together. From Tables 1 and 2 it is clear that the cooling part is much easier to replicate accurately than the heating one. Almost any combination gives deviations less than 20% in the cooling for any of the stars, while only about half of them minimise χ_H below 20% for the F, G, and K and only a few for the M star. Table 3 shows how many cases we find with both χ_C and χ_H smaller than certain thresholds indicated individually for each of the stars. For the M2V star the choice of four τ bins is obviously insufficient. There are only a few combinations (less than 1%) that lead to χ_C < 10% and χ_Η < 20% and there is not a single one that reduces χ_C and χ_Η below 5 and 20% simultaneously.

To understand how the deviation measures vary with the location of the three separators, an illustrative example for the G2V star is shown in Fig. 10. The figure shows the dependence of χ_C and χ_Η on the location of the middle separator τ₁ when the other two separators are fixed. To visualise the results, in each column of the figure the deepest separator is fixed at one of six values (log τ₀ = 0.5, 0, −0.5, −1.2, −1.9, −2.9). Each of the curves correspond to a different fixed value of the highest separator τ₂. The position of the highest separator is indicated by the colour, the black curve corresponding to the smallest value of τ₂ and the red curve to the largest one.

The figure shows that χ_C and χ_H vary significantly with the location of τ₁ if the location of the deepest separator τ₀ is deeper than some critical depth (left-hand columns of Fig. 10). That critical location of τ₀ is located deeper for χ_C than for χ_H (log τ₀ ≃ −0.5 for χ_C and log τ₀ ≃ −1.2 for χ_H). For χ_C especially, if τ₀ ≃ −0.5 the curves flatten out, meaning that the results are essentially insensitive to the location of the other two nodes.

If τ₀ is deeper than these critical values, the deviation measures are insensitive to the location of the other two separators only if τ₂ is relatively close to τ₀ (the red curve in the left-hand column of Fig. 10). When τ₀ 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 τ₀ ≥ −0.5, the shape of the curves for the χ_C and χ_H in Fig. 10 becomes more different. The cooling deviation shows one minimum around log τ₁ ≈ −2.5, while the heating deviation has two minima: one for the locations of τ₁ close to τ₀ and another one that moves higher up with increasing height of τ₀. When log τ₀ = −0.7, the minimum is at log τ₁ ≈ −5, and when log τ₀ ≥ −0.5, the minimum is at log τ₁ ≈ −2.5. Since the location of one of the minima of χ_H coincides with the minimum of χ_C, both can be minimised using the same set of separators even when τ₀ is as deep as log τ₀ = 0.5.

When the location of the separator τ₀ is set above the critical values, the sensitivity of both χ_C and χ_H to the location of τ₁ and τ₂ mostly disappears (see the sequence from the left to the right panel of Fig. 10). This happens for a higher up location of τ₀ for χ_H than for χ_C^.

Our analysis suggests that, in the case of the G2V star, there are two strategies for placing the separators so that χ_C and χ_H remain below 4 and 15% respectively. The first strategy is to set log τ₀ around the continuum formation height, [−0.5, 0] and log τ₁ ∈ [−3, −2.5]. The second one is to place log τ₀ so that it includes most of the weak and intermediate photospheric lines (log τ₀ ∈ [−2.4, −1.9]) and log τ₁ even higher in [−5, −3.5]. In both cases the location of log τ₂ can be anywhere above log τ₁. The former case (see 2nd and 3rd column) contains the combinations that produce the minimal χ values in our experiment, although the curves are not entirely flat in such a way that the optimal solution is sensitive to the model atmosphere and to the ODF data. Alternatively, the latter case still guarantees small deviations, but it appears to be more robust as χ_C and χ_H essentially depend only on the location of the deepest bin. Since our model atmosphere does not have a chromospheric temperature rise, however, the strategy with all separators positioned relatively high may be limited to that type of models.

For the spectral types F3V and K0V the results are similar to those for G2V. For τ₀ located deeper than a certain critical value, χ_C and χ_H vary with τ₁ similarly to the G star, showing two minima for χ_H, of which one coincides with the single minimum of χ_C. There is a little variation of the deviation curves with τ₂ except when it is close to the τ₀. For the F3V star, the optimal combinations that produce χ_C ≤ 8% and χ_H ≤ 15% are (a) log τ₀ ∈ [−0.2, 0.5] and log τ₁ ∈ [−2.5, −2] and (b) log τ₀ ∈ [−2.2, −0.7] and log τ₁ ∈ [−5.5, −4.7], with log τ₂ again anywhere above log τ₁ .

For the K0V star, we identified three variants that produce χ_C ≤ 5% and χ_H ≤ 15%: (a) log τ₀ ∈ [0, 0.3] and log τ₁ ∈ [−2, −1], (b) log τ₀ ∈ [−0.5, −0.2] and log τ ∈ [−4, −3], and (c)log τ₀ ∈ [−1.9, −1.7] and log τ₁ ∈ [−5, −3], with log τ₂ < log τ₁.

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 τ₀ there is only one minimum for the χ_H varying with τ1 and it does not coincide anymore with the minimum for χ_C, but with its local maximum (see two left columns in Fig. 11). Secondly, once the τ₀ separator is sufficiently high in the atmosphere so that the deviation measures are insensitive to changes in τ₁, χ_C flattens out as for the other stars but its value quickly increases with the height of τ₀ (see sequence from left to right in the figure), from χ_C ≃ 6% when log τ₀ = 0.2 to χ_C ≃ 15.5% when log τ₀ = −0.9. Higher up, the value of χ_C remains larger than 12%, while χ_H also starts to increase once it becomes flat (after τ₀ = −1.2). Therefore, for the M2V star it is not possible to minimise both deviation measures simultaneously with four bins. The best combination of parameters that we found gives χ_C ≤ 9% and χ_H ≤ 27% for log τ₀ ∈ [0.2, 0.5] and log τ₁ ∈ [−1.1, −0.6], with log τ₂ < log τ₁.

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 τ = {−1.9, −5, −7} for the F3V star; log τ = {0, −2.5, −5} for the G2V; log τ = {−0.2, −3, −5} for the K0V; log τ = {0.3, −0.8, −2.5} for the M2V.

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

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

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

5.2 Number of τ-bins

To understand the dependence of the Q rate on the number of separators, we computed Q for the four stars for different number of separators (2, 3, 5, 10, 20, 30, 40, 50, and 100)⁷ that are uniformly distributed in the relevant range for each of the stars: log τ^ref ∈ [−9, 0.5] for F3V, [−6.5, 0.5] for K0V, and G2V and [−4.5, 0.5] for M2V. Again, the selected optical depth ranges ensure that all the stellar models are properly represented (see the vertical axis of Fig. 8).

In Fig. 12 χ_C and χ_H are shown as a function of the number of bins. For F-, G-, and K - type stars, there is very little variation in χ_C The trend is slightly different for each of the stars individually. While χ_C slowly decreases for the K star with an increasing number of bins, it increases for G2V, and it first increases and then decreases for the F3V – type star. In all three cases the value of χ_C saturates at around 20 bins. The M star shows a significant improvement in χ_C when the number of bins increases from 4 to 21. For more bins it saturates at ≈5%. Although there is more variation in χ_Η than in χ_C, the overall trend is similar and, again, the results saturate for more than ≈20 bins. Vögler et al. (2004) tested the difference between the Q rate computed with the binned opacity and with the ODF for an increasing number of bins (up to 60). Similarly to our results, they found that the solutions for Q converge towards a limiting value that is different from Q_ODF.

Wray et al. (2006) investigated the performance of the OB and ODF methods in the context of 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 , ϰ-sorting instead of τ-sorting), they also found that the results of the radiation computation quickly converge with increasing number of bins (reaching saturation at around 10 bins). However, note that their conclusions about the ODF solution are not directly translatable to the radiative transfer in the stellar atmospheres owing to different definitions of ODF.

Therefore, the saturation values of χ_C and χ_H are the minimal deviations that cannot be eliminated by increasing the number of the τ-separators. Depending on the stellar type, the minimal χ_C values span between 2 and 7%, and the minimal χ_H span between ≈10 and ≈30%. These minimal deviations are the smallest for the K and the largest for the F star. However, it is possible to find certain combinations of the separators that minimise the deviations particularly well. For example, the χ_H value for the F3V star is reduced to only a few per cent when three separators are set to log τ = {−1.8, −4.3, −6.6}. This is consistent with locations that produce simultaneous minimum χ_C and χ_H deviations, as discussed in Sect. 5.1. The deviations with the optimal combinations found using four bins for each of the stars in Sect. 5.1 are shown as coloured plus markers in Fig. 12. These deviations are either similar to or smaller than the saturation values when the bins are distributed uniformly. We have also computed the Q rates using the Rosseland opacity mean applied to the whole spectrum (Eq. (C.1)). The χ differences for this grey solution are shown as coloured circles in Fig. 12. Their values are larger or similar to the saturation ones, but always larger than the solutions with four bins from Sect. 5.1.

We experimented with several other strategies for distributing separators automatically. They all confirm the conclusion that careful distribution of the separators is more important than their number.

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

5.3 {τ, λ} -binning

The OB method using only τ-separators does not reproduce the exact solution in the limit of infinite number of bins because the contributions from different parts of the spectrum, each with its own height variation, are mixed in any given bin (Vögler et al. 2004, see discussion in Appendix C). One way to improve the approximation is to preserve some of the wavelength dependence in the binning procedure by grouping the opacity data points with respect to their spectral regions prior to the τ-binning. Therefore, λ-separators between these regions are introduced in addition to the r-separators. Such bidimensional distribution of the opacities is non-uniform in (τ, λ) parameter space: in every wavelength interval the τ-separators may be specified independently (both their number and location) to account for the height distribution of the opacities in that interval. Ludwig (1992) introduced the {τ, λ}-binning, showing that splitting the continuum at the Balmer jump improved significantly the results for Vega, due to the systematically different behaviour of Q for wavelengths shortwards and longwards of the jump. Vögler (2003) explored the possibility of splitting the least opaque bin into two subgroups for long and short wavelengths and concluded that the improvement is probably not worth the extra computational effort. Collet et al. (2018) also attempted to improve the OB procedure by partly accounting for the lost wavelength dependence. Magic (2014) treated the UV and the visible + IR opacities separately. Each group is sorted with respect to its own set of τ-separators. In addition, the least opaque bin of the visible + IR opacities is further split with a set λ-separators.

We computed the Q values using a {τ, λ}-binning approach similar to Magic (2014, see his Fig. 2.6) for different total number of bins, N = 6, 7, 8, 9, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, and 100. The distribution of the bins for the G2V star for N = 6, 7, 8 and 9 is shown in Fig. 13. For any N, the spectrum is first divided into the short and long wavelength parts by placing one λ-separator at λ₀. The value of λ₀ 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 τ-bins, and the long-wavelength ones into two. The least opaque τ-bin of the long-wavelength opacities is then split by one or more additional λ-separators. For N = 6 the λ-separator is placed at 1600 nm. For N > 6, the λ-separators of the least opaque τ-bin are placed equidistantly in log λ in the range [λ₀, 9570] nm. The τ-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 τ-λ plane (cf. Figs. 8 and 9).

We computed the Q rate for every combination of {τ, λ}-bins. The corresponding deviations χ_C and χ_H are shown in Fig. 14. The behaviour of χ_C with the number of bins is generally similar to that found in Fig. 12. For four stars the χ_C deviation decreases with increasing number of bins until it saturates after around 10–20 bins. The saturation values of χ_C are lower than the ones found in Sect. 5.2 (cf. Table 4). For less than 15 bins (Δ log λ > 0.1) χ_C remains nearly constant for G2V and K0V, but increases significantly for F3V and M2V. The results for M2V show two peaks for N = 7 and 10, and a minimum between them for N = 8.

This strong variation in χ_C is a good illustration of how sensitive the Q rate is on how the ODF data points are distributed between the bins. Figure 13 shows the distribution of the bins with N = 7 and N = 9 (2nd and 4th panel). Bins 1 to 5 are identical in both cases. Bins 6 and 7 for N = 7 are split into two bins each for N = 9. Bin 7 (for N = 7) includes most of the visible spectrum, it coincides with the central part of the Planck function and thus it dominates the radiative losses. The rate Q computed from that bin is compared to the total Q in the 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 χ_C, this leads to a drop from ≈14% for N = 7 to ≈3% for N = 9 (top panel of Fig. 14). The finer the distribution of the visible opacities in wavelength is, the more stable are the results leading to the saturation of χ_C.

For the M2V star the χ_H saturation value is about the same as in Fig. 12, but for the F, G, and K stars χ_H saturates at significantly larger values owing to an excess of heating relative to the ODF solution.

In another experiment, we tested a variant of the previous method were bins are distributed in a checkerboard pattern in the (τ, λ) plane. The number of τ-separators is varied between 1 and 16 and the number of λ-separators is fixed at 5 (λ = 381, 726, 1383, 2636, 5023 nm). The τ-separators for all the stars are distributed equidistantly at log τ^ref ∈ [−3, 0.5]. The minimum number of bins in this experiment is 12, the maximum 102. The results are shown in Fig. 16. For a small number of bins χ_H is larger than the deviation measure we obtain with the optimised solutions from Sect. 5.1. With 12 bins the values χ_H lie between 20 and 30%, similar to those in Fig. 12 with 11 bins distributed only in τ. With finer sampling in τ, the value of χ_H is quickly reduced below 10% and saturated at that level for all four stars.

From this analysis it follows that introduction of the λ separators may reduce the minimum deviation of the OB in particular cases when it is combined with carefully selected separators in τ. However, this improvement is relatively small in comparison with the basic method where sorting is done only by optical depths (cf. Figs. 16 and 12). Table 4 summarises the deviation measures χ_C and χ_H computed in the experiments described in this section. The values shown for the experiments described in Sects. 5.2 and 5.3 correspond to the minimum number of bins at which the deviation measures approach their saturation values. It is clear that the 4-bin approach with carefully tuned locations of the τ-separators (Sect. 5.1) performs similarly or better than the other approaches for all four stars. The only exception is the M star, where taking into account the wavelength distribution simultaneously with distribution in τ for all λ intervals is an obvious advantage (Fig. 14). The strategy proposed by Magic et al. (2013) reproduces the cooling component better than other strategies, but it produces large errors in the heating component.

Finally, in Fig. 17 Q and Q/ρ are compared for the four stars and five opacity approaches (ODF, grey opacity, and three OB setups). The deeper part of the atmosphere where mainly cooling occurs is shown in the two left columns; the higher layers where the heating occurs are shown in the two right columns. On one hand, the cooling can be fitted either using Q or Q/ρ, since the shape of Q is practically not modified by the relatively slowly changing density in 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/ρ show similar information regarding the quality of the computed Q^OB in respect to Q^ODF. Some differences are expected when choosing Q/ρ, but the general behaviours described in this work should be the same. If one is interested in layers closer to the chromosphere and beyond, Q/ρ would be the quantity to optimise.

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

Fig. 15 Comparison of the radiative energy exchange rate Q from bin 7 for N = 7 and bins 8 and 9 for N = 9 corresponding to the cases shown in Fig. 14 for the M2V star. Left panel: Q7 dominates the total QOB = Σl Ql for N = 7. Right panel: comparison of Q8 and Q9 for N = 9 against Q7 for N = 7.

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

6 Conclusions

This work is based on detailed monochromatic opacity tables computed with the SYNSPEC (Hubeny & Lanz 2011) code using complete 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 τ bins (Sect. 5.2) and simultaneously binning in optical depth and wavelength (Sect. 5.3). None of the strategies tested converge to the ODF solution in the limit of a large number of bins since the binning mixes parts of the spectrum with different height profiles of the opacity (Appendix C). It is possible to find combinations of small number of bins that significantly reduce the deviations. However, these combinations are model-dependent and with relatively narrow intervals in which the bins can be located. Strategies that include preservation of the wavelength dependence through combined τ- and λ-binning may improve the deviation measures, but with a significant increase in computing overheads. For example, even in the most extreme case of the M star, to reduce the already small χ values by a factor of three (see the bottom row of Table 4), requires an increase in the number of bins, and thus the computing cost by a factor of four or five. In most applications such accuracy is not needed as other approximations made in the RTE solution and the MHD modeling contribute more to the total error budget.

The results presented in this study are limited to the solar chemical composition. In the stars with different metallicities both the atmospheric structure and the relative importance of the various opacity contributors will change. Therefore, our results should not be blindly applied to these stars. Nevertheless, the trends identified in our analysis may be used as guidelines for tuning the details of the OB procedure for the stars with non-solar metallicities.

An open question left in this work is the importance of the measured deviations χ_C and χ_H in the final structure of the atmosphere in the context of simulations. To answer this question, one may start by extending the study for the Sun in Sect. 3.4 of Vögler et al. (2004) for the case of other stellar spectral types in time dependent simulations.

It should also be emphasised that the optimal binning approximation selected based on the 1D models will not necessarily remain optimal when implemented in 3D simulations. However, the tuning of the binning strategy is computationally too expensive to be routinely done in 3D. In the 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.

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

Acknowledgements

This work was supported by the European Research Council through the Consolidator Grant ERC–2017–CoG–771310–PI2FA and by Spanish Ministry of Science through the grant PID2021–127487NB–I00. A.P.G. acknowledges support from the Agencia Estatal de Investigación (AEI) of the Ministerio de Ciencia, Innovación y Universidades (MCIU) and the European Social Fund (ESF) under grant with reference PRE2018–086567. C.A.P. is thankful for funding from the Spanish goverment through grants AYA2014–56359–P, AYA2017–86389–P and PID2020–117493GB–100. We thank the referee for the detailed and careful reading of the manuscript and for asking the questions that lead to Appendix C and Fig. 17. We thank Hans Günter Ludwig for reading the manuscript, the interesting discussions about the OB method and for suggesting the analysis in Fig. C.4. We are thankful to Alexander Shapiro for reading the manuscript and giving his feedback. We are thankful to Terry Mahoney for the language editing. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

Appendix A Content of the SYNSPEC monochromatic opacity

In SYNSPEC each of the transitions taken into account to compute the opacity can be specified by the user. Only explicit ions (in our computations, H, H⁺, H⁻, He, He⁺, He⁺⁺ C, C⁺, C⁺⁺, N, N⁺, N⁺⁺, O, O⁺, O⁺⁺, Na, Na⁺, Na⁺⁺, Mg, Mg⁺, Mg⁺⁺, Al, Al⁺, Al⁺⁺, Si, Si⁺, Si⁺⁺, Ca, Ca⁺, Ca⁺⁺, Fe, Fe⁺ and Fe⁺⁺) are considered in the continuum opacities solved in LTE (see sections 2.4 and 2.5 of Hubeny & Lanz 2017a as well as Hubeny et al. 2021 for more details). We use the data included in the synple repository (version 1.2), as described in the following paragraphs. The contributions involving molecules and any of the Collision-Induced Absorptions are only taken into account below certain temperature, in this work chosen to be 8000K (i.e. when molecules are considered in the EOS).

The 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).

• from Gingerich (1969).

• 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 TOP-base⁸ and the Iron Project ⁹:

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

• : 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 hydro-genically 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 H₂ Rayleigh scattering from Dalgarno (1962).

• Thomson scattering: ϰ = 6.65 × 10⁻²⁵n_e/ρ (cm²g⁻¹)

• Collision-Induced Absorption of H₂-H₂ (Borysow et al. 2001), H₂-He (Jørgensen et al. 2000), H₂ -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 ¹⁰, updated by data from the National Institute of Standards and Technology (NIST) when available. SYNSPEC also uses data from EXOMOL¹¹ for some of the molecular transitions. In the case of the molecules, around a total of 2 × 10⁷ transitions for H₂, CH, NH, OH, NaH, MgH, SiH, CaH, CrH, FeH, C₂, CN, CO, MgO, AlO, SiO, CaO, VO, H₂O and TiO are included (see Fig. 5). In the case of the atoms, SYNSPEC includes a line list with around 2 × 10⁶ transitions for:

• Neutral H, As, Se, Rb, Sb, Te, Cs, Pt, Au, Tl, and Bi.

• Neutral and first ion of He, Li, Ga, Ge, Sr, Y, Tc, Pd, Ag, Cd, In, Sn, Ba, La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Tm, Yb, Lu, Hf, Ta, W, Re, Os, Ir, Hg, Pb, Th, and U.

• Neutral and first two ions of Be, Nb, Rh, Zr, and Ru.

• Neutral and first three ions of B, C, and Mo.

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

• Neutral and first seven ions of Zn.

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

• Neutral and first nine ions of Fe, Ni, and Cu.

Appendix B Details of the comparison of the opacity distribution function from SYNSPEC and ATLAS

Figures in this appendix show examples of typical differences between the ODFs from SYNSPEC and ATLAS in the opacity and in the Q rates for the four stellar models (see Sect. 4.2).

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

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

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

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

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

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

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

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

Appendix C Saturation of the opacity binning method and comparison against grey opacity

The OB method is not expected to converge to the ODF solution for a large number of τ-bins (Vögler et al. 2004). To explain why this is the case, we use the checkerboard distribution of the bins introduced in Fig. C.1, similar to that used in Sect. 5.3, but now with four bins in τ and six in λ. The number of non-empty τ-bins changes depending on the λ-bin. For example, in the case of the F3V star there are four τ-bins for λ < 300 nm, while there are two τ-bins for λ ∈ [1000, 2000] nm (Fig. C.1). For each star, in Fig. C.2 we compare the binned opacity of pairs of wavelength ranges selected to better illustrate this behaviour. We compare λ-ranges that have the same number and location of non-empty τ-bins. For each star and any of the τ-bins, the opacity from one of the two λ-ranges vary with height in a different way than the opacity from the other. Even for arbitrarily increased number of bins, binning only in τ mixes wavelength regions which opacities have different height profiles (this is illustrated in figure 4.13 from Vögler 2003, too). Although this effect may be reduced, it is present even when splitting in λ, since each contributor to the opacity may have its own variation with height.

Figure C.3 shows how the opacity computed with the OB method converges to the Rosseland opacity in each τ-bin, (dotted black lines), in the limit of large enough optical depths. We show only the case for the F3V star and the range for λ < 300 nm to avoid overcrowding the plot. However, the conclusions are valid for any λ-range and star. Another example of this convergence can be found in figure 4.12 from Vögler (2003) for the case of the Sun. Figure C.3 also shows the Rosseland opacity mean for the whole spectrum in black dashed line. The harmonic mean used to compute this averaged opacity,

(C.1)

gives more weight to the continuum than to the line cores, making that the Rosseland mean is close to the opacity of the least opaque bins.

Finally, Figure C.4 shows that the combination of the opacity bins in the Rosseland fashion,

(C.2)

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

(C.3)

thus, the combination of the bins following Eq. C.2 gives

(C.4)

This equation is similar to Eq. C.1, but with the formalism of the ODF to compute the integrals and with the additional sums over the bins. It can be written in a more general way as

(C.5)

where λ(l) corresponds to the wavelengths that enter the lth bin. Binning in τ, λ, or both always is equivalent to group only separating in λ, since each opacity point corresponds to a different wavelength¹². This property as well as understanding the integrals of the opacity in terms of Riemann or Lebesque integration explains the result shown in Fig. C.4. Similarly, the Planck mean should be recovered for higher layers.

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

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

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

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

References

All Tables

All Figures

	Fig. 1 Temperature stratification of the F3V (green), G2V (yellow). K0V (orange), and M2V (dark red) 1D hydrostatic models used in this work. The same colour coding is used throughout the paper to label results for the four stellar types. The grid used to solve the RTE is marked by dots.
In the text

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

In the text

Fig. 3 Monochromatic radiative energy exchange rate Qv (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 |Qv| to max |Qv|. Blue tones in the colour-map correspond to negative Qv, i.e. cooling in the atmosphere, and red tones correspond to positive Qv, i.e. heating.

In the text

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

In the text

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

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

In the text

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

In the text

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

In the text

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

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

In the text

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

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

In the text

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

In the text

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

In the text

	Fig. 15 Comparison of the radiative energy exchange rate Q from bin 7 for N = 7 and bins 8 and 9 for N = 9 corresponding to the cases shown in Fig. 14 for the M2V star. Left panel: Q7 dominates the total QOB = Σl Ql for N = 7. Right panel: comparison of Q8 and Q9 for N = 9 against Q7 for N = 7.
In the text

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

In the text

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

In the text

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

In the text

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

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

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

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

In the text

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

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

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

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

In the text

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

In the text

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

In the text

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