Korg: a modern 1D LTE spectral synthesis package

We present Korg, a new package for 1D LTE (local thermodynamic equilibrium) spectral synthesis of FGK stars, which computes theoretical spectra from the near-ultraviolet to the near-infrared, and implements both plane-parallel and spherical radiative transfer. We outline the inputs and internals of Korg, and compare synthetic spectra from Korg, MOOG, Turbospectrum, and SME. The disagreements between Korg and the other codes are no larger than those between the other codes, although disagreement between codes is substantial. We examine the case of a C$_2$ band in detail, finding that uncertainties on physical inputs to spectral synthesis account for a significant fraction of the disagreement. Korg is 1-100 times faster than other codes in typical use, compatible with automatic differentiation libraries, and easily extensible, making it ideal for statistical inference and parameter estimation applied to large data sets. Documentation and installation instructions are available at https://ajwheeler.github.io/Korg.jl/stable/.


INTRODUCTION
Improvements in instrumentation have yielded exponential growth in the amount of spectral data to analyse. Creating pipelines that can keep up with analysis is nontrivial. There are several extant codes for 1D LTE spectral synthesis, including Turbospectrum (Plez et al. 1993;Plez 2012;Gerber et al. 2022), Moog (Sneden 1973;Sneden et al. 2012), SYNTHE (Kurucz 1993;Sbordone et al. 2004), SME, (Valenti & Piskunov 1996, 2012Piskunov & Valenti 2017;Wehrhahn et al. 2022), SPECTRUM (Gray & Corbally 1994), and SYNSPEC (Hubeny & Lanz 2011Hubeny et al. 2021). While they have enabled a huge volume of research, these codes can be difficult to use for the uninitiated, and require input and output through custom file formats, impeding integration into analysis code. Here we present Korg, a new 1D LTE synthesis package, written in Julia and suitable for easy integration with scripts and use in an interactive environment. As the first such new code in more than two decades, Korg benefits from numerical libraries not available at the time earlier packages were authored, principally modern automatic differentiation packages and optimization libraries.
The two fundamental assumptions made by Korg are that the stellar atmosphere is hydrostatic and 1D, and in LTE (local thermodynamic equilibrium). Eliminating the former two assumptions means calculating atmospheric structure using the full hydrodynamic equations (e.g. Freytag et al. 2012;Magic et al. 2015;Schultz et al. 2022), or the magnetohydrodynamic equations (e.g. Vögler et al. 2005), if the internal magnetic field is strong.
Fortunately, corrections to 1D LTE level populations (e.g. Amarsi et al. 2020Amarsi et al. , 2022, equivalent widths, or abundances (e.g. Lind et al. 2011;Amarsi et al. 2015;Bergemann et al. 2012;Osorio & Barklem 2016) can be calculated from NLTE simulations. These can be applied to LTE results or codes to produce approximate NLTE spectra at relatively little computational cost. Additionally, biases from the assumption of LTE will roughly cancel for similar stars, yielding differential abundance estimates with high precision.
The performance of spectral synthesis codes is most important when fitting observational data. Because synthesis must be embedded in an inference loop, the analysis of of a single spectrum may trigger tens or hundreds . Vertical offsets have been applied, but the relative scaling is preserved, which is why some derivative spectra appear flat here (but show fluctuations at a sub-percent level).
of syntheses. Even when many parameters may be estimated by interpolating over (or otherwise comparing to) a precomputed grid of spectra (e.g. Recio-Blanco et al. 2006;Smiljanic et al. 2014;Holtzman et al. 2018;Boeche et al. 2021;Buder et al. 2021), individual abundances are best done with targeted syntheses of individual lines. Furthermore, generating a grid from which to interpolate can be computationally expensive. Inference and optimization can also be sped up by fast and accurate derivatives of the function being sampled or minimized, most easily produced via automatic differentiation. Korg is designed to be compatible with automatic differentiation packages (e.g. ForwardDiff), which can provide derivative spectra in roughly the same amount of time required for a single synthesis (as discussed in Section 4; see Figure 1).

DESCRIPTION OF CODE
To synthesize a spectrum, Korg calculates the number density of each species (e.g. H I, C II, CO) at each layer in the atmosphere (Section 2.2), then computes the absorption coefficient at each wavelength and atmospheric layer due to continuum (Section 2.3) and line absorption (Section 2.4). Given the total absorption coefficient at each wavelength and atmospheric layer, it then solves the radiative transfer equation to produce the flux at the stellar surface (Section 2.5).

Inputs
Korg takes as inputs, a model atmosphere, a linelist, and abundances for each element in A(X) form, 1 assumed to be constant throughout the atmosphere. Korg includes functions for parsing .modformat MARCS (Gustafsson et al. 2008) model atmospheres, and line lists in the format supplied by VALD (Piskunov et al. 1995;Kupka et al. 1999Kupka et al. , 2000Ryabchikova et al. 2015;Pakhomov et al. 2019) or Kurucz (Kurucz 2011), or those accepted by Moog. When parsing VALD linelists, Korg will automatically apply corrections to unscaled log(gf ) values for transitions with isotope information using the isotopic abundance values supplied by Berglund & Wieser (2011) for linelists which 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 linelists to vacuum (Birch & Downs 1994). Users may also construct line lists or model atmospheres using custom code and pass them directly into Korg.

Chemical equilibrium
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 non-LTE (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 2 (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, 3, the correction is smaller than 10 −3 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).
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 which 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.

Continuum absorption
Korg computes contributions to the continuum absorption coefficient from a number of sources, listed in Table 1. Continuum absorption mechanisms involving 2 https://physics.nist.gov/PhysRefData/ASD/levels form.html  an atom and an electron can be classified as bound-free (bf) or free-free (ff), depending on whether the elec- tron is initially bound to the atom. 3 In addition to bound-free and free-free absorption for H − , H I, He II, and 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 4 , when available, TOPBase 5 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 wasn't 3 Note that these are named as though the electron were bound, even for free-free interactions. For example, "H I bf" refers to the ionization of neutral hydrogen, and "H I ff" refers to absorption by the interaction of free electrons and free protons. 4 https://norad.astronomy.osu.edu/#AtomicDataTbl1 5 https://cds.u-strasbg.fr/topbase/topbase.html present. We considered all species with bound-free interaction included in Gustafsson et al. (2008), but included only those which contribute to the total continuum absorption at above the 10 −3 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. 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 linelist in Korg. The agreement between the continua calculated by Korg and other codes is good (see Section 3).

Line absorption
The contribution of each line to the total absorption coefficient, α line , is where the wavelength-integrated cross-section, σ is given by 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, λ 0 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 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 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 linelist, Korg approximates it with 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 de-exitation dominates the transition's energy uncertainty, and that the upper level's degeneracy is unity. The values of γ stark and γ vdW at 10 5 K, γ stark 0 and γ vdW 0 , are provided by the linelist, then scaled by their temperature dependence according to semiclassical impact theory (e.g. Hubeny & Mihalas 2014 ch. 8.3) to the per-particle broadening parameters: 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 linelist, 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, r 2 , is approximated by 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 set to the value at which a pure Gaussian or pure Lorentzian profile takes on a value of α crit , whichever is larger. By default, Korg truncates profiles at α crit = 10 −3 α 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 (Stehlé & Hutcheon 1999), which are pre-convolved with Doppler profiles. For H α , H β , and H γ , where selfbroadening is important, we add to φ a Voigt profile using the p-d approximation for self-broadening from Barklem et al. (2000b).

Radiative transfer
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 in a plane-parallel atmosphere, and in a spherical atmosphere, where I λ is the intensity of radiation at wavelength λ, S λ the source function, α λ the absorption coefficient, z the negative depth into the atmosphere, r, the distance to the center of the star, and µ 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 it's 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, Here, I top λ (µ) 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 λ = (R/d) 2 F λ , 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, where E 2 is the second order exponential integral, τ λ is the optical depth (dτ λ = α λ dz) and τ λ 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, F λ , Korg calculates the intensity, I top λ at a discrete grid of µ values, by integrating along rays from the lowest atmospheric layer to the top atmospheric layer for many of surface µ values. This is valid since the source function is isotropic under LTE. Rays which 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). The agree at the sub-percent level in flux, and complication with the Bezier scheme lead us to make the lower-order scheme the default. Appendix A contains numerical details and an extended discussion of the calculation of transfer integrals.

COMPARISON TO OTHER CODES
To test Korg's correctness, we compare its spectra to those produced by Moog, Turbospectrum, and SME for four sets of the parameters (T eff , log g, and metallicity). The first set of parameters is solar, and we have chosen the others to be similar to three well-studied stars: the red giant branch star Arcturus, HD122563 (an extremely metal-poor red giant), and HD49933 (an F-type main sequence star), but shifted slightly to the nearest existing MARCS model atmosphere to avoid the need to interpolate the atmospheres. For all stars, the same solar abundance pattern was assumed, which is not problematic since we are comparing synthetic spectra to each other and not to observational data. For the first four regions, we used an "extract all" linelist from VALD, which includes all known lines within the wavelength bounds. The near-infrared spectrum was synthesized using the APOGEE DR16 linelist (Smith et al. 2021). The VALD linelists are pre-adjusted for isotopic abundances, which eliminates a possible source of disagreement. For the APOGEE linelist, Turbospectrum's runtime isotopic adjustment may not be consistent with Korg's or with the pre-adjustment performed for SME or Moog, but agreement is generally good nevertheless. As Moog has a limit of 2500 lines which can contribute at a given wavelength, we culled the weakest TiO lines from each list, reducing the total number of transitions passed to Moog to be at most 10,000. This is expected to have negligible impact on the spectrum. Since Korg (see Section 2.3) and Moog (excepting the fork presented in Sobeck et al. 2011) do not support true scattering, we set the PURE-LTE flag to true in babsma.par, setting Turbospectrum to turn off this functionality. Scattering can't be turned off in SME which presumably explains some of the disagreement between it and other codes. Disagreement between the codes is substantial, with fluxes disagreeing at the ten percent level in many cases, as established by Blanco-Cuaresma (2019). Agreement is strongest in the infrared, where continuum opacities are simpler. We discuss some of the discrepancies in more detail here.

Continuum absorption in the violet and ultraviolet
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 defined structure seen the Korg continuum is due to resonances in the metal 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-effects of line blanketing).

Balmer series
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 shal-low features, in contrast to Turbospectrum (hydrogen lines are omitted from the Moog linelist in this region). They are located at 3668.76Å, 3667.17Å, 3665.75 A, 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 6 (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.
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 Stehlé & Hutcheon (1999) selfbroadening 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 constrast, 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 modelling of the H α core requires techniques beyond 1D LTE (e.g. Barklem 2007;Amarsi et al. 2018), and that the observational data doesn't match the prediction of any of the codes, we consider this disagreement permissible.

C 2 band
While tracking down the causes of disagreement between codes for all wavelengths and stars isn't feasible, it is instructive to focus on one as a "case study". In the sun, Korg predicts deeper lines from roughly 5160 A -5165Å than the other codes. These are due to absorption by C 2 , and the disagreement arises from the varying molecular equilibrium constants, adopted for the species by the codes. Most of the difference in K C2 values comes from different values for the dissociation energy, D 0 , of C 2 . Molecular dissociation energies are one of the most dominant sources of systematic uncertainty in spectral synthesis (see discussion in Barklem & Collet 2016). To summarize the treatments of K C2 in each code: • Korg uses the data from Barklem & Collet (2016): D 0 = 6.371 eV.
• Based on logging information, Turbospectrum appears to use the polynomial expansions from Tsuji (1973): estimated D 0 = 6.07 eV.
• The Moog molecular equilibrium constants comes from Kurucz 7 . We were able to extract the polynomial approximation for the C 2 molecular equilibrium constant from the source code: D 0 = 6.21 eV.
• For C 2 , the SME equilibrium constant comes from Sauval & Tatum (1984): D 0 = 6.297 eV. Figure 10 shows the effects of these choices about K C2 . They account for nearly all of the disagreement in the synthesized spectra. While the numerical scheme used to represent the temperature-dependence of K C2 plays some role, the value of D 0 adopted is far more important. The differing D 0 values result in number densities of C 2 differing by up to a factor of 2 at the photosphere, and by orders of magnitude at the top of the atmosphere.

BENCHMARKS
The first panel in Figure 11 shows the time taken by each code to produce the spectra in the first four wavelength regions in Section 3. These are relatively small regions with linelists including all transitions in VALD, including many which have essentially no effect on the spectrum. The second panel shows the time to compute  Barklem & Collet (2016) data adjusted for the D0 used by another code. Uncertainty in D0 drives most of the disagreement in the amount of C2 present at the photosphere and most of the disagreement in the synthesized spectra at these wavelengths.
a spectrum from 15000Å -15500Å using the APOGEE DR16 linelist. This provides a more realistic example of synthesizing a spectrum as part of the analysis for a large survey. All test are single-core, and run on the same machine with an AMD Epyc 7702P.
Korg only needs to load the linelist and model atmosphere into memory once, so repeat syntheses which use the same inputs (as would be the case when varying individual abundances) avoid that step. For these comparisons, loading the linelist (10 5 -10 6 transitions) takes 1 -4 seconds, and loading the model atmosphere takes roughly 0.05 ms. As nearly all cases requiring performance involve repeated synthesizing with the same (or nearly the same) linelist, we did not optimize input parsing for speed and do not include this time in Korg benchmarks. We note that for the "small wavelength regions", the Moog times are technically lower limits, since the linelists provided to Moog were reduced in size by an order of magnitude or more (though the synthesized spectra are essentially unchanged). Figure 11 has separate markers for Turbospectrum with and without hydrogen lines, as we found them to slow down synthesis by a very large factor. In wavelength ranges where hydrogen lines are not present or unimportant, they can be omitted, and Turbospectrum executes much more quickly. The comparison to SME is included for completeness, but is not entirely fair since SME includes the effects of scattering, which is numerically expensive. Scattering is turned off for Turbospectrum, to make the comparison as fair as possible.
In the top panel of Figure 12 we plot the time required by Korg to compute the gradient of the solar spectrum, ∂F ∂A(X) , with regard to a varying numbers of element abundances, N . (See Figure 1 for a subset of the derivative spectra.) For this demonstration we used the 15,000 -15,500Å range and the APOGEE DR16 linelist. A dashed horizontal line marks the time required to synthesize the spectrum, just under one second. Calculating the N -dimensional gradient spectrum takes roughly 2 + 0.15N seconds, meaning that the marginal time required for each derivative is an order of magnitude smaller than the time required to obtain it via finite differences. As Julia's automatic differentiation ecosystem improves, Korg may benefit from further speedups to the calculation of gradient spectra. Calculating derivative spectra with respect to atmospheric parameters, e.g. T eff , log g, or v mic , requires code to interpolate between model atmospheres, which we plan to add to Korg.

CONCLUSIONS AND FUTURE DEVELOPMENT
We have presented a new code, Korg, for 1D LTE spectral synthesis, i.e. computing stellar spectra given a model atmosphere, linelist, and element abundances. Korg is both fast faster than all other codes, and autodifferentiable, which yields further speedups when synthesis is embedded in an optimization loop. The code is publicly available at https://github.com/ajwheeler/ Korg.jl and installable via the Julia package manager. Detailed documentation, usage examples, and installation instructions are available at https://ajwheeler. github.io/Korg.jl/stable.
In comparing Korg to other codes, we have highlighted that the level of disagreement between them (and thus systematic error) is substantial. Much of this can be attributed to uncertain physical parameters (e.g. the C 2 dissociation energy, discussed in Section 3.3) and

Sol
Arcturus HD122563 HD49933  Figure 11. The time taken by each code for the syntheses in Section 3. top: average compute time for the optical wavelength regions, which used very dense linelists. bottom: compute time to synthesize spectra from 15,000 -15,500Å with the APOGEE DR16 linelist, a more "realistic" use case. Note that the wavelength regions synthesized were larger than those plotted in Figures 5 -8. SME includes a true scattering treatment which can't be turned off, increasing its execution times. models of continuum absorption processes, but varying numerical methods also play a role (e.g. the discussion in appendix A). Unfortunately, similar problems extend to model atmospheres, which require the same physics used for synthesis, and linelists, which include many poorlyconstrained atomic parameters.
We have shown in this work that Korg produces results for FGK stars that are similar to other codes. We recommend its use for these spectral classes. We plan to soon address the factors limiting Korg's applicability outside this regime, most importantly its lack of polyatomic molecules (present in cool stars) and true scattering (important in hot stars). We also plan to add support for NLTE departure coefficients.
There are also limitations affecting ease-of-use that we plan to ad-0 1 2 3 4 5 6 7 8 9 10 gradient dimensions (element abundances) dress. As mentioned in Section 4, calculating derivative spectra with regard to atmospheric parameters like log g or T eff , requires code to interpolate model atmospheres which is written in Julia. Relatedly, while inferring stellar abundances and parameters with Korg is not overly burdensome, the user must still apply the line spread function to synthesized spectra and calculated goodness-of-fit themselves, processes we would like to automate.
There are some limitations which apply to all spectroscopic synthesis codes. In the blue and ultraviolet, the poorly-constrained continuum results in very uncertain predictions, manifesting in the poor agreement found between codes. Many of the chemical equilibrium constants necessary for spectroscopic synthesis are also poorly-constrained, as discussed Barklem & Collet (2016). There are doubtless many features (like the one discussed in Section 3.3) where chemical equilibrium constants are the primary driver of uncertainty. Uncertain atomic parameters in linelists pose a similar problem in all regimes.
Our goal is that Korg will be useful for both inference of stellar parameters and abundances for large survey data, e.g. Gaia  Zhao et al. 2012), GALAH (De Silva et al. 2015), and for boutique analyses of individual spectra. We aim to make Korg fast and flexible enough to enable better survey pipelines and novel analyses, such as the propagation of error from synthesis inputs to synthesized spectra, or the joint inference of line parameters with observational data. In addition, we hope that by making Korg as easy to use as possible, more researchers will find it worthwhile to synthesize spectra when the need arises. 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 solution (mτ − b) exp(−τ )dτ = − exp(−τ ) (b + m(τ + 1)) .
Calculating the optical depth, τ , requires numerically integrating where α 5 is the absorption coefficient at the reference wavelength, 5000Å, calculated by Korg, and α 5 is the absorption coefficient at the reference wavelength supplied by the model atmosphere. Their ratio is the correction factor to α λ which enforces greater consistency between the model atmosphere and spectral synthesis. Korg computes τ λ by evaluating the equivalent integral where τ 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 ln τ 5 , and the integrand of A4 is more nearly linear than that of A3.

A.2. Comparison to quadratic Bézier scheme
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 non-monotonic 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 A.2 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 sub-percent 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.

S)
Bezier spline Figure 13. The difference in the rectified spectrum obtained for the Sun using the quadratic Bézier scheme and Korg's default scheme. The purple and brown lines show the result obtained using modified Bézier interpolation, and cubic spline interpolation of α λ , respectively. (This is separate from the interpolation of the source function, which uses Bézier curves in both cases.) The difference between all three methods is small compared to the difference in spectra obtained from different codes, but nevertheless larger than is ideal.