Korg: A Modern 1D LTE Spectral Synthesis Package

Korg takes as inputs, a model atmosphere, a line list, and abundances for each element in A(X) form, ⁶ assumed to be constant throughout the atmosphere. Korg includes functions for parsing .mod-format MARCS (Gustafsson et al. 2008) model atmospheres, and line lists in the format supplied by VALD (Piskunov et al. 1995; Kupka et al. 1999, 2000; Ryabchikova et al. 2015; Pakhomov et al. 2019) or Kurucz (Kurucz 2011) or those accepted by Moog. When parsing VALD line lists, Korg will automatically apply corrections to unscaled $\mathrm{log}({gf})$ values for transitions with isotope information using the isotopic abundance values supplied by Berglund & Wieser (2011) for line lists that are not already adjusted for isotopic abundances. Korg detects and uses packed Anstee, Barklem, and O'Mara (ABO; Anstee & O'Mara 1991; Barklem et al. 2000a) if they are provided (as VALD optionally does). It uses vacuum wavelengths internally and will automatically convert air-wavelength VALD line lists to vacuum (Birch & Downs 1994). Users may also construct line lists or model atmospheres using custom code and pass them directly into Korg.

Korg uses the assumption of LTE, i.e., that in a sufficiently small region, baryonic matter is described by thermal distributions, and radiation is only slightly out of detailed balance. The source function (the ratio of the per-volume emission and absorption coefficients) is dominated by collisions of baryonic matter and is Planckian. While NLTE calculations are important for producing the most unbiased possible model spectra, they are prohibitively slow and not yet suitable for applications that require computing many spectra over large wavelength ranges.

In order to compute the absorption coefficient, α_λ , at each layer of the atmosphere, Korg must solve for the number density of each species. Treating the number density of each neutral atomic species as the free parameters, Korg uses the NLsolve library (Mogensen et al. 2020) to solve the system of Saha and molecular equilibrium equations with the temperature, total number density, and number density of free electrons set by the model atmosphere, and the total abundance of each element set by the user. By default, Korg uses molecular equilibrium constants from Barklem & Collet (2016; the ionization energies are originally from Haynes 2010), but alternatives can be passed by the user. These are defined only up to 10,000 K, so we treat the molecular partition functions as constant above this temperature. This is unlikely to be problematic as few molecules are present above this threshold. The default atomic partition functions are calculated using energy levels from NIST ⁷ (Kramida & Ralchenko 2021).

In dense environments, upper energy levels are perturbed or dissolved. This can be crudely accounted for by truncating the terms of the partition function (e.g., Hubeny & Mihalas 2014, Section 4.1) or with the "probability occupation formalism" developed in Hummer & Mihalas (1988) and generalized by Hubeny et al. (1994). This effect is most important for hydrogen, where it strongly impacts the partition function starting above 10,000 K. Based on energy level truncation, the partition functions of lithium and vanadium are affected in the same regime at the 5% and 3% level, respectively, but other elements are unchanged at the 1% level. For FGK stars, these effects are most important for higher-order hydrogen transitions. Figure 2 shows the occupation probability correction factor, w, for several energy levels of hydrogen in the solar atmosphere. For n = 1, 2, and 3, the correction is smaller than 10⁻³ everywhere. At present, we do not include these effects in Korg because they do not make a large difference for FGK stars and because of disagreement with observational data (see Section 3.2).

Zoom In Zoom Out Reset image size

By default, Korg solves the equilibrium equations taking into account all elements up to uranium (as neutral, singly ionized, and doubly ionized species), as well as 247 diatomic molecules. As they exist in extremely small quantities in LTE atmospheres, Korg neglects species that are more than doubly ionized. It also neglects ionic and triatomic molecules, which are present in cool stars, although we plan to address this in a future version.

Korg computes contributions to the continuum absorption coefficient from a number of sources, listed in Table 1. Continuum absorption mechanisms involving an atom and an electron can be classified as bound–free (bf) or free–free (ff), depending on whether the electron is initially bound to the atom. ⁸ In addition to bound–free and free–free absorption for H⁻, H i, He ii, and ${{\rm{H}}}_{2}^{+}$ (as well as He⁻ and He), Korg includes treatments of bound–free and free–free interactions with metals. Free–-free interactions are treated with the hydrogenic approximation using Gaunt factors from van Hoof et al. (2014), with corrections from Peach (1970) for neutral He, C, Si, and Mg. Bound–-free interactions are treated with cross sections from NORAD, ⁹ when available, and TOPBase ¹⁰ (Seaton et al. 1994) otherwise. Following Gustafsson et al. (2008), we shifted the theoretical energy levels to the empirical values from NIST, leaving out levels for which an empirical counterpart was not present. We considered all species with bound–free interaction included in Gustafsson et al. (2008), but included only those that contribute to the total continuum absorption above the 10⁻³ level anywhere in any of the atmospheres used in Section 3. Figure 3 shows the contribution of various mechanisms to the continuum absorption coefficient at the solar photosphere. At present, Korg treats all scattering as absorption, which is correct in the LTE regime. In the future we plan to support quasi-LTE radiative transfer with isotropic scattering, which will yield more accurate spectra when Rayleigh scattering dominates or is a significant source of opacity.

Zoom In Zoom Out Reset image size

Table 1. Sources of Continuum Absorption in Korg

Absorption Source	Wavelength Bounds	Temperature Bounds	Reference
H− bf	≥1250 Å		McLaughlin et al. (2017)
H− ff	2604–113, 918 Å	2520–10,080 K	Bell & Berrington (1987)
He− ff	5063–151,878 Å	2520–10,080 K	John (1994)
H i bf			Karzas & Latter (1961) via Gingerich (1969) via Kurucz (1970)
H i ff	100–106 Å	100–106 K	van Hoof et al. (2014)
He ii bf			Karzas & Latter (1961) via Gingerich (1969) via Kurucz (1970)
${{\rm{H}}}_{2}^{+}$ bf and ff	700–20,000 Å	3150–25,200 K	Stancil (1994)
He ii ff	100–106 Å	100– 106 K	van Hoof et al. (2014) with departure coefficients from Peach (1970)
metals ff	100– 106 Å	100– 106 K	van Hoof et al. (2014) with departure coefficients from Peach (1970) for C i, Si i, and Mg i
C i bf	500–30000 Å	100– 105 K	Nahar & Pradhan (1991)
Na i bf	500–30000 Å	100– 105 K	Seaton et al. (1994)
Mg i bf	500–30000 Å	100– 105 K	Seaton et al. (1994)
Al i bf	500–30000 Å	100– 105 K	Seaton et al. (1994)
Si i bf	500–30000 Å	100– 105 K	Nahar & Pradhan (1993)
S i bf	500–30000 Å	100– 105 K	Seaton et al. (1994)
Ca i bf	500–30000 Å	100– 105 K	Seaton et al. (1994)
Fe i bf	500–26088 Å	100– 105 K	Bautista (1997)
H Rayleigh	≥1300 Å		Lee (2005) via Colgan et al. (2016)
He Rayleigh	≥1300 Å		Dalgarno & Kingston (1960), Dalgarno (1962), Schwerdtfeger (2006) via Colgan et al. (2016)
H2 Rayleigh	≥1300 Å		Dalgarno & Williams (1962)
electron scattering			Thomson (1912)

Note. When bounds in temperature or wavelength are missing, the absorption coefficient is defined for all positive values.

Download table as:  ASCIITypeset image

Figure 4 shows the ratio of the continuum absorption according to Korg and the one used by MARCS through the solar atmosphere. Agreement is good, except in the violet and ultraviolet, where Korg uses more recent bound–free metal absorption coefficients with more structure in wavelength, which shows up as vertical bands in the figure. Compounding this, the MARCS values also include the effects on line blanketing, which is handled via the line list in Korg. The agreement between the continua calculated by Korg and other codes is good (see Section 3).

Zoom In Zoom Out Reset image size

The contribution of each line to the total absorption coefficient, α_line, is

$$\begin{eqnarray}&&{\alpha }_{\lambda }^{\mathrm{line}}=\sigma n\displaystyle \frac{\exp (-\beta {E}_{\mathrm{up}})-\exp (-\beta {E}_{\mathrm{lo}})}{U(T)}\phi (\lambda ),\end{eqnarray} \tag{ 1 }$$

where the wavelength-integrated cross section, σ, is given by

$$\begin{eqnarray}&&\sigma ={g}_{\mathrm{lo}}f\displaystyle \frac{\pi {e}^{2}{\lambda }_{0}^{2}}{{m}_{e}{c}^{2}}.\end{eqnarray} \tag{ 2 }$$

Here, n is the number density of the line species, E_up and E_lo are the upper and lower energy levels of the transition, β = 1/kT is the thermodynamic beta, ϕ is the normalized line profile, U is the species partition function, T is the temperature, f is the oscillator strength, g_lo is the degeneracy of the lower level, λ₀ is the line center, and m_e is the electron mass.

For all lines of all species besides hydrogen (including autoionizing lines), ϕ is approximated with a Voigt profile using the numerical approximation from Hunger (1956). The width of the Gaussian component is

$$\begin{eqnarray}&&{\sigma }_{D}={\lambda }_{0}\sqrt{\displaystyle \frac{2{kT}}{m}+{\xi }^{2}},\end{eqnarray} \tag{ 3 }$$

where m is the species mass, and ξ is the microturbulent velocity, a fudge factor used to account for convective Doppler shifts. The width of the Lorentz component (in frequency rather than wavelength units) is given by

$$\begin{eqnarray}&&{\rm{\Gamma }}={{\rm{\Gamma }}}_{\mathrm{rad}}+{\gamma }_{\mathrm{stark}}{n}_{e}+{\gamma }_{\mathrm{vdW}}{n}_{{\rm{H}}\ {\rm\small{}}{\rm{i}}},\end{eqnarray} \tag{ 4 }$$

where Γ_rad is the radiative broadening parameter, γ_stark the per-electron Stark broadening parameter, γ_vdW the per-neutral-hydrogen van der Waals broadening parameter, and n_e and n_{H I} are the number densities of free electrons and neutral hydrogen, respectively. We neglect pressure broadening of molecular lines, setting γ_stark and γ_vdW to zero.

When Γ_rad is not supplied in the line list, Korg approximates it with

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{rad}}=\displaystyle \frac{8{\pi }^{2}{e}^{2}}{{mc}{\lambda }^{2}}{gf},\end{eqnarray} \tag{ 5 }$$

where m is the mass of the atom or molecule, g is the degeneracy of the lower level of the transition, and f is the transition oscillator strength. This can be obtained by assuming that spontaneous deexitation dominates the transition's energy uncertainty, and that the upper level's degeneracy is unity.

The values of γ_stark and γ_vdW at 10⁵ K, ${\gamma }_{0}^{\mathrm{stark}}$ and ${\gamma }_{0}^{\mathrm{vdW}}$, are provided by the line list, then scaled by their temperature dependence according to semiclassical impact theory (e.g., Hubeny & Mihalas 2014, chap. 8.3) to the per-particle broadening parameters:

$$\begin{eqnarray}&&{\gamma }_{\mathrm{stark}}={\gamma }_{0}^{\mathrm{stark}}{\left(\displaystyle \frac{T}{{10}^{5}\,{\rm{K}}}\right)}^{\tfrac{1}{6}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\gamma }_{\mathrm{vdW}}={\gamma }_{0}^{\mathrm{vdW}}{\left(\displaystyle \frac{T}{{10}^{5}\,{\rm{K}}}\right)}^{\tfrac{3}{10}}.\end{eqnarray} \tag{ 7 }$$

If provided, ABO parameters, which describe a more nuanced temperature dependence in broadening by neutral hydrogen, will be used to calculate γ_vdW instead. When a Stark broadening parameter is not provided in the line list, the approximation from Cowley (1971) is used. When a van der Waals broadening parameter is provided, a form of the the Unsöld approximation (Unsold 1955; Warner 1967) is used in which the angular momentum quantum number is neglected and the mean square radius, $\overline{{r}^{2}}$, is approximated by

$$\begin{eqnarray}&&\overline{{r}^{2}}=\displaystyle \frac{5}{2}{\left(\displaystyle \frac{{{n}_{\mathrm{eff}}}^{2}}{Z}\right)}^{2},\end{eqnarray} \tag{ 8 }$$

where n_eff is the effective principal quantum number and Z is the atomic number. Pressure broadening is neglected for autoionizing lines with no provided parameters.

The absorption coefficient for each line (except those of hydrogen) is calculated over a dynamically determined wavelength range. The maximum detuning (wavelength difference from the line center) is calculated by taking the quadrature sum of the detunings at which a pure Gaussian and pure Lorentzian profile reach the value of α_crit. Korg truncates profiles at α_crit = 10⁻³ α_cntm, where α_cntm is the local continuum absorption coefficient. This ratio, α_crit/α_cntm, can be set by the user.

The broadening of hydrogen lines is treated separately. We use the tabulated Stark broadening profiles from Stehle & Hutcheon (1999), which are pre-convolved with Doppler profiles. For H_α , H_β , and H_γ , where self-broadening is important, we add to ϕ a Voigt profile using the p-d approximation for self-broadening from Barklem et al. (2000b).

Given the absorption coefficient, α_λ , at each wavelength and atmospheric layer, the final step is to solve the radiative transfer equation. The radiative transfer equation is

$$\begin{eqnarray}&&\displaystyle \frac{\mu }{{\alpha }_{\lambda }}\displaystyle \frac{{\mathrm{dI}}_{\lambda }}{\mathrm{dz}}={S}_{\lambda }-{I}_{\lambda },\end{eqnarray} \tag{ 9 }$$

in a plane-parallel atmosphere, and

$$\begin{eqnarray}&&\displaystyle \frac{\mu }{{\alpha }_{\lambda }}\displaystyle \frac{\partial {I}_{\lambda }}{\partial r}+\displaystyle \frac{1-{\mu }^{2}}{{\alpha }_{\lambda }r}\displaystyle \frac{\partial {I}_{\lambda }}{\partial \mu }={S}_{\lambda }-{I}_{\lambda },\end{eqnarray} \tag{ 10 }$$

in a spherical atmosphere, where I_λ is the intensity of radiation at wavelength λ, S_λ is the source function, α_λ is the absorption coefficient, z is the negative depth into the atmosphere, r is the distance to the center of the star, and μ is the cosine of the angle between the r/z axis and the line of sight. When the thickness of the atmosphere is small relative to the stellar radius, curvature can be neglected and a plane-parallel atmosphere is a good approximation; otherwise, sphericity must be taken into account. By default, Korg does its radiative transfer calculations in the same geometry as the model atmosphere.

What Korg actually returns is the disk-averaged intensity, i.e., the astrophysical flux,

$$\begin{eqnarray}&&{{ \mathcal F }}_{\lambda }=2\pi {\int }_{0}^{1}\mu {I}_{\lambda }^{\mathrm{top}}(\mu )\,{\rm{d}}\mu .\end{eqnarray} \tag{ 11 }$$

Here, ${I}_{\lambda }^{\mathrm{top}}(\mu )$ stands for I_λ (z = 0, μ) in the plane-parallel case and I_λ (r = R, μ) in the spherical case, where R is the radius of the outermost atmospheric layer. The total flux from the star is then ${F}_{\lambda }={(R/d)}^{2}{{ \mathcal F }}_{\lambda }$, where d is the distance to the star. Since Korg assumes that the stellar atmosphere is in LTE, S_λ is the blackbody spectrum.

In the plane-parallel case,

$$\begin{eqnarray}&&{{ \mathcal F }}_{\lambda }=2\pi {\int }_{{\tau }_{\lambda }^{\prime} }^{0}{S}_{\lambda }({\tau }_{\lambda }){E}_{2}({\tau }_{\lambda }){\rm{d}}{\tau }_{\lambda },\end{eqnarray} \tag{ 12 }$$

where E₂ is the second-order exponential integral, τ_λ is the optical depth (dτ_λ = α_λ dz), and ${\tau }_{\lambda }^{\prime}$ is the depth at the bottom of the atmosphere. The spherical case admits no further analytic simplifications. When using a spherical model atmosphere, to obtain the astrophysical flux, ${{ \mathcal F }}_{\lambda }$, Korg calculates the intensity, ${I}_{\lambda }^{\mathrm{top}}$, at a discrete grid of μ values, by integrating along rays from the lowest atmospheric layer to the top atmospheric layer for many surface μ values. This is valid since the source function is isotropic under LTE. Rays that do not intersect the lowest atmospheric layer are cast from the far side of the star. The integral over μ (Equation (11)) is performed using Gauss–Legendre quadrature.

Korg has two radiative transfer calculation schemes: a lower-order default and the quadratic Bezier scheme from de la Cruz Rodríguez & Piskunov (2013). They agree at the subpercent level in flux, and complications with the Bezier scheme lead us to make the lower-order scheme the default. The Appendix contains numerical details and an extended discussion of the calculation of transfer integrals.

Near the Balmer jump and blueward, agreement between the synthetic continua is generally poor. This is the primary cause for the disagreement between the codes in the first two wavelength regions, which are at the shortest wavelengths. The disagreement at the Balmer jump between Turbospectrum and the other codes arises from the fact that it dissolves the outer hydrogen orbitals according to the Hummer & Mihalas (1988) formalism when calculating the H i bf absorption coefficient. The defined structure seen in the Korg continuum is due to resonances in the bf cross sections, not lines. Unfortunately, it is the case that reliable spectral synthesis in the blue and ultraviolet remains out of reach, as metal opacities are theoretically very challenging to predict at the energy resolution required (not to mention the difficult-to-model effects of line blanketing).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are eight Balmer series lines in the 3660–3670 Å (vacuum) wavelength region (upper levels n = 26 to n = 33), five of which are included in Korg (those up to n = 30). Typically, these have a minimal impact on the observed spectrum, but in the metal-poor giant, Korg predicts that these lines are visible as broad shallow features, in contrast to Turbospectrum (hydrogen lines are omitted from the Moog line list in this region). They are located at 3668.76, 3667.17, 3665.75, 3664.48, and 3663.33 Å (vacuum), and are more clearly visible in the residuals panel. This is because Turbospectrum uses the occupation probability formalism discussed in Section 2.2, which eliminates the relevant hydrogen orbitals throughout most of the stellar atmospheres. In a high-resolution ultraviolet spectrum of HD122563 from the Ultraviolet and Visual Echelle Spectrograph ¹¹ (Dekker et al. 2000), a star with similar parameters, the lines are present. For this reason, we are not confident that the occupation probability formalism offers a significant advantage over unmodified partition functions, and we do not use it in Korg. The top left panel in Figure 9 shows part of this wavelength region in more detail, along with the observational data. SME predicts the strongest absorption by these lines in the metal-poor giant, but it also predicts strong absorption for the other stars, where it is not present in observations. The bottom left panel of Figure 9 demonstrates this for the Sun, using the Wallace et al. (2011) solar atlas.

Zoom In Zoom Out Reset image size

There is relatively good agreement between Korg, Turbospectrum, and SME in the H_α wings (Moog lacks special treatment of hydrogen lines and generally disagrees). Turbospectrum and SME both use forms of HLINOP (Barklem & Piskunov 2003), so strong agreement is to be expected (though they presumably have differing treatments of the equation of state). Korg uses roughly the same treatment as HLINOP, with Stark broadening from Stehle & Hutcheon (1999) and self-broadening via the Voigt approximation from Barklem & Piskunov (2003), but Figure 9 indicates that there are differences in implementation. For the metal-poor giant hot dwarf (Figures 7 and 8), Turbospectrum and SME agree with each other almost exactly and disagree with Korg at the ∼1% level in the H_α wings. In contrast, Korg and Turbospectrum agree more closely with each other than they do with SME for the Sun. In the H_α core, disagreement is similarly minor, except for in the metal-poor giant. The right panel in Figure 9 shows this in detail, alongside the observed spectrum of HD122563, a star with similar parameters, from GALAH (Buder et al. 2021). Given that accurate modeling of the H_α core requires techniques beyond 1D LTE (e.g., Barklem 2007; Amarsi et al. 2018), and that the observational data does not match the prediction of any of the codes, we consider this disagreement permissible.

While tracking down the causes of disagreement between codes for all wavelengths and stars is not feasible, it is instructive to focus on one as a "case study." In the Sun, Korg predicts deeper lines from roughly 5160–5165 Å than the other codes. These are due to absorption by C₂, and the disagreement arises from the varying molecular equilibrium constants,

$$\begin{eqnarray}&&{K}_{{{\rm{C}}}_{2}}=\displaystyle \frac{{P}_{{\rm{C}}}^{2}}{{P}_{{{\rm{C}}}_{2}}},\end{eqnarray} \tag{ 13 }$$

adopted for the species by the codes. Most of the difference in ${K}_{{{\rm{C}}}_{2}}$ values comes from different values for the dissociation energy, D₀, of C₂. Molecular dissociation energies are one of the most dominant sources of systematic uncertainty in spectral synthesis (see the discussion in Barklem & Collet 2016). To summarize the treatments of ${K}_{{{\rm{C}}}_{2}}$ in each code:

• 1.  
  Korg uses the data from Barklem & Collet (2016): D₀ = 6.371 eV.
• 2.  
  Based on logging information, Turbospectrum appears to use the polynomial expansions from Tsuji (1973): estimated D₀ = 6.07 eV.
• 3.  
  The Moog molecular equilibrium constants comes from Kurucz. ¹² We were able to extract the polynomial approximation for the C₂ molecular equilibrium constant from the source code: D₀ = 6.21 eV.
• 4.  
  For C₂, the SME equilibrium constant comes from Sauval & Tatum (1984): D₀ = 6.297 eV.

Figure 10 shows the effects of these choices about ${K}_{{{\rm{C}}}_{2}}$. They account for nearly all of the disagreement in the synthesized spectra. While the numerical scheme used to represent the temperature dependence of ${K}_{{{\rm{C}}}_{2}}$ plays some role, the value of D₀ adopted is far more important. The differing D₀ values result in number densities of C₂ differing by up to a factor of 2 at the photosphere and by orders of magnitude at the top of the atmosphere.

Zoom In Zoom Out Reset image size

To solve the equation of radiative transfer along a ray, Korg approximates the source function, S (λ subscripts dropped for brevity), with linear interpolation over optical depth, τ. Between each adjacent pair of atmospheric layers, we have S(τ) ≈ m τ + b, so the (indefinite) transfer integral has the solution

$$\begin{eqnarray}&&\int (m\tau -b)\exp (-\tau ){\rm{d}}\tau =-\exp (-\tau )\left(b+m(\tau +1)\right).\end{eqnarray} \tag{ A1 }$$

When evaluating Equation (12), Korg uses the same approximation,

$$\begin{eqnarray}\begin{array}{rcl}\int (m\tau -b){E}_{2}(\tau ){\rm{d}}\tau & = & \displaystyle \frac{1}{6}\left[\tau {E}_{2}(\tau )(3b+2m\tau )\right.\\ & & \left.-\,\exp (\tau )(3b+2m(\tau +1))\right].\end{array}\end{eqnarray} \tag{ A2 }$$

Calculating the optical depth, τ, requires numerically integrating

$$\begin{eqnarray}&&{\tau }_{\lambda }(s)={\int }_{0}^{s}{\alpha }_{\lambda }\displaystyle \frac{\alpha {{\prime} }_{5}}{{\alpha }_{5}}\mathrm{ds},\end{eqnarray} \tag{ A3 }$$

where α₅ is the absorption coefficient at the reference wavelength, 5000 Å, calculated by Korg, and $\alpha {{\prime} }_{5}$ is the absorption coefficient at the reference wavelength supplied by the model atmosphere. Their ratio is the correction factor to α_λ that enforces greater consistency between the model atmosphere and spectral synthesis. Korg computes τ_λ by evaluating the equivalent integral

$$\begin{eqnarray}&&{\tau }_{\lambda }(s)={\int }_{0}^{\mathrm{ln}\tau {{\prime} }_{5}(s)}\tau {{\prime} }_{5}\displaystyle \frac{{\alpha }_{\lambda }}{{\alpha }_{5}}{\rm{d}}(\mathrm{ln}\tau {{\prime} }_{5}),\end{eqnarray} \tag{ A4 }$$

where $\tau {{\prime} }_{5}$ is the optical depth at the reference wavelength per the model atmosphere. This integral is numerically preferable, since atmospheric layers are spaced uniformly in $\mathrm{ln}\tau {{\prime} }_{5}$, and the integrand of (A4) is more nearly linear than that of (A3).

Korg has a secondary radiative transfer scheme based on de la Cruz Rodríguez & Piskunov (2013), which approximates the source function with monotonic quadratic Bézier interpolation. They suggest that τ_λ might be computed using the same scheme, but this causes numerical problems when the absorption coefficient is nonmonotonic in depth. We modified the interpolant by requiring the Bézier control point (Equation (10) in de la Cruz Rodríguez & Piskunov 2013) to be within a reasonable range, which eliminates the issue, but is physically unmotivated. This is the implementation available in Korg. As an alternative, we tried computing τ_λ by interpolating α_λ with a cubic spline. Unfortunately, the relatively sparse sampling of model atmospheres means that this choice has an impact on the synthesized spectrum similar to the difference between the default method and the de la Cruz Rodríguez & Piskunov (2013) method. Figure A1 shows the difference in a portion of the solar spectrum resulting from using Bézier interpolation of the source function compared to the default implementation. The difference in synthesized spectra is at the subpercent level, significantly smaller than the level of disagreement between codes. While we plan to adopt a higher-order scheme as the default in the future, we provide this scheme as a secondary option until we better understand the impact of the choices involved.

Zoom In Zoom Out Reset image size