Korg: Fitting, Model Atmosphere Interpolation, and Brackett Lines

We describe several updates to Korg, a package for 1D LTE spectral synthesis of FGKM stars. Built-in functions to fit observed spectra via synthesis or equivalent widths make it easy to take advantage of Korg's automatic differentiation. Comparison to a past analysis of 18 Sco shows that we obtain significantly reduced line-to-line abundance scatter with Korg. Fitting and synthesis are facilitated by a rigorously tested model atmosphere interpolation method, which introduces negligible error to synthesized spectra for stars with T eff ≳ 4000 K. For cooler stars, atmosphere interpolation is complicated by the presence of molecules, though we demonstrate an adequate method for cool dwarfs. The chemical equilibrium solver has been extended to include polyatomic and charged molecules, extending Korg's regime of applicability to M stars. We also discuss a common oversight regarding the synthesis of hydrogen lines in the infrared, and show that Korg's Brackett line profiles are a much closer match to observations than others available. Documentation, installation instructions, and tutorials are available at https://github.com/ajwheeler/Korg.jl.


INTRODUCTION
Korg, first presented in Wheeler et al. (2023), is a stellar spectral synthesis code for computing spectra of FGKM stars.It it been written with both flexibility and performance in mind, and is faster than similar codes by a factor of 1-100.It's central assumptions are of 1D geometry and local thermodynamic equilibrium, which, while not state-of-the art theoretically, are well understood and significantly speed up and simplify the calculation of spectra.Korg is written in pure Julia, which makes it easy to use in scripts or an interactive programming environment from either Julia or Python.It also supports automatic differentiation, which can speed up calculation involving derivatives by an order of magnitude or more.
This paper describes several updates to Korg which aim to make it useful to researchers in more contexts (Wheeler et al. 2023 provides an overview of the code, here we focus on changes only).Section 2 discusses Corresponding author: Adam Wheeler wheeler.883@osu.edunew routines for automatic fitting of observed data via synthesis or equivalent widths, as well as the employed model atmosphere interpolation scheme.This requires the interpolation of model atmospheres, which is not often given the attention it requires.Section 3 discusses several model improvements, most notably a more sophisticated treatment of molecules (section 3.1), which allows for synthesis of M star spectra.Section 4 summarizes our conclusions.

FITTING SPECTRA VIA DIRECT COMPARISON OR EQUIVALENT WIDTHS
Korg now includes functions for inferring parameters and abundances via direct fitting of the observed spectrum (synthesis) or via equivalent widths, find best fit params and ews to abundances.

Fitting via synthesis
The find best fit params function is designed to make it as easy as possible to take advantage of Korg's automatic differentiation capabilities, even when called from Python.It takes as input an observed rectified spectrum with known error, a linelist and initial guesses for each parameter to be fit.Fitting is performed with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (see, e.g.Nocedal & Wright 2006) as implemented in the Optim.jlpackage (Mogensen et al. 2020).The fit is performed within wavelength windows specified by the user, with minimal computational overhead imposed by non-contiguous windows.Any subset of the following parameters can be fit, with the remainder fixed to default or user-specified values: • T eff : the effective temperature.
• log g: the surface gravity.
• m H: the metallicity of the star, in dex.(defaults to 0.0).
• alpha H: the α-element abundance of the star, in dex.(defaults to the value of m H.) Here, the αelements are interpreted as those with even atomic numbers from C to Ti.
• Individual element abundances, e.g.'Na', specify the solar-relative ([X/H]) abundance of that element (defaults to the value of the metallicity) .
• vsini: the projected rotational velocity of the star, in km/s (defaults to 0.0).
• epsilon: the linear limb-darkening coefficient (defaults to 0.6).This is used only for applying rotational broadening.
By default, convolution with the line-spread function (LSF) is automatically handled via construction of a sparse matrix, but user-specified LSFs can be substituted if desired.

Fitting via equivalent widths
The ews to abundances function computes abundances of individual atomic transitions, given some measured equivalent widths and a model atmosphere.This functionality is inspired by the widely used ABFIND routine in MOOG (Sneden 1973), but differs in its approach.While MOOG computes a fictitious curve of growth for each species, Korg synthesizes each line in wavelength segments and integrates the rectified profile to compute a model equivalent width.Nearby lines are segmented into groups to ensure wavelength segments do not overlap, and that absorption from nearby lines do not contribute.Chemical equilibrium and continuum opacity calculations are not duplicated within groups of wavelength segments to avoid redundant computation.
Individual line abundances are computed by first synthesizing each line at an initial abundance 1 specified by the user, [A(X)] 0 , to obtain an equivalent width, W 0 .Lines are assumed to be on the linear part of the curve of grown, so then, where W is the observed equivalent width and A(X) is the resulting abundance.While this procedure could be sped up by using a fictitious curve of growth for each line, synthesizing each line once ensures the most accurate calculation possible.Figure 1 shows an example use of ews to abundances functionality.Equivalent widths are measured and reported by Meléndez et al. (2014) for the Sun and 18 Sco, a solar-like twin.Korg was used to interpolate a model atmosphere for the Sun (T eff = 5777 K, log g = 4.44, [Fe/H] = 0.0) and for 18 Sco, using the parameters reported by Meléndez et al. (2014): T eff = 5823 K, log g = 4.45, [Fe/H] = 0.054.Because Korg's model atmosphere grid (discussed in section 2.3) does not include v mic , both model atmospheres were interpolated with the default value for dwarfs, v mic = 1 km s −1 .Subsequent syntheses of each profile used v mic = 1 km s −1 for the sun, and v mic = 1.02 km s −1 for 18 Sco (the same as in Meléndez et al. 2014).The differential (18 Sco -Sun) abundances for neutral (circle) and singly ionized (square) transitions are shown as a function of excitation potential and reduced equivalent width.
While the same stellar parameters and equivalent widths are used, the scatter in Korg line abundances is about half (σ = 0.0056 compared to σ = 0.010) that reported in Meléndez et al. (2014).The weakest lines show the largest abundance differences between 18 Sco and the Sun, where even the few significant digits used to record the equivalent width (e.g., rounding from 12.15 m Å to 12.2 m Å) can translate to a systematic shift of 0.01 dex for that line abundance.Figure 2 compares the line-to-line scatter for all abundances measured using equivalent widths in Meléndez et al. (2014), demonstrating that the improved precision is not limited to iron lines.Excluding the spectral synthesis code used, the other notable differences between these analyses are the grid of model atmospheres (Meléndez et al. 2014 used Kurucz models), how they are interpolated (section 2.3), and synthesizing lines instead of using a fictitious curve of growth.Thus, Korg's ew to abundances functionality is suitable for classical spectroscopic deter-  Figure 2. The line-to-line scatter in the abundance of each element measured using equivalent widths from Meléndez et al. (2014).Comparison between Korg reanalysis and original abundances.Computing abundances with Korg results in significantly smaller scatter for most elements.mination of stellar parameters, including high-precision differential-abundance research.

Model atmosphere interpolation
In order to synthesize a spectrum for a star with arbitrary parameters, or to fit a observed spectrum, the input model atmosphere must be either constructed ad hoc or interpolated to precise parameter values.Korg now includes functionality to interpolate the Sloan Digital Sky Survey (SDSS) model atmosphere grid2 , an expanded version of the MARCS (Gustafsson et al. 2008) grid documented in Mészáros et al. (2012).The grid contains 579,150 model atmospheres, which vary over five parameters: T eff (2500 to 8000 K), log g (-0.5 to 5.5), [metals/H] (-2.5 to 1.0), [α/Fe] (-1 to 1), and [C/Fe] (-1 to 1).Atmospheres for some parameters (mostly in unphysical regimes, but also the tip of the red-giant branch) did not converge.Mészáros et al. (2012) filled these in with radial basis function interpolation.The atmospheres in the grid are divided into two groups.Those with log g > 3 are in plane-parallel geometry and were computed with a microturbulent velocity of 1 km s −1 .Those with log g < 3 are in spherical geometry (with an assumed mass of 1M ⊙ ) and were computed with a microturbulent velocity of 2 km s −1 .
When Korg takes model atmosphere as input, it uses five quantities defined at each layer of the atmosphere: temperature, number density, electron number density (used only as a starting guess, see section 3.1), physical depth, and the optical depth at 5000 Å.We experimented with various transforms of the parameters and found that better accuracy is obtained when number density and electron number density are log trans-formed, and when geometric depth is transformed with sinh −1 .Korg interpolates the transformed quantities as a function of the parameters over which the grid varies using simple multilinear interpolation (provided by the Interpolations.jlpackage3 ).See the discussion below regarding cubic interpolation of model atmospheres for cool stars.
We use a simple procedure to quantify errors introduced by model atmosphere interpolation.First we construct a new, temporary grid by interpolating halfway between nodes in all dimensions, then we interpolate back to the the original grid points.Differences between the rectified spectra produced from a doubly interpolated atmosphere and the corresponding original atmosphere provides an upper bound on the error introduced by interpolation.
Figure 3 shows the largest absolute error introduced to syntheses in the visible and infrared by the interpolation of model atmospheres in all five parameters.Interpolations errors are slightly lower in the infrared than in the visible because, as noted in Westendorp Plaza et al. ( 2023), traditional interpolation methods perform the worst for the highest and lowest atmospheric layers, and lines in the infrared are less sensitive to the upper atmosphere.For most parts of the Kiel diagram, the error introduced by interpolation is negligible, but it becomes important for cool stars.Because we are comparing rectified spectra, the largest errors are in the line cores, though for cool line-blanketed stars, the interpolation error manifests as a near-constant offset.In contrast with Mészáros & Allende Prieto (2013), we find that there is no benefit to performing interpolation at the spectral level, rather than the atmospheric level (which is more costly unless spectra are precomputed).
While model atmosphere interpolation works well for most stars, it is problematic for those with T eff ≲ 4000 K.As an experiment, we developed an alternative interpolation method (not yet available in Korg) tailored for cool dwarfs.We found that τ 5000 , the optical depth at 5000 Å, is particularly difficult to accurately interpolate in the cool-dwarf regime.When using the default radiative transfer scheme, Korg uses τ 5000 to anchor its calculations (see Wheeler et al. 2023).The method tailored for cool dwarfs performs interpolation on atmospheres that have been resampled (via cubic interpolation) so that atmospheric depth is now indexed in terms of a shared set of fixed τ 5000 .In other words, we no longer interpolate τ 5000 between grid points.For cool dwarfs, resampling the atmospheres had negligible effect on synthesized spectra at grid points and significantly reduced the interpolation error (estimated as described above).Additionally, we interpolated between the resampled model atmospheres using cubic (rather than linear) interpolation.The results of this procedure are shown in Figure 4.Both resampling and cubic interpolation improve the interpolation error for cool dwarfs, reducing the error to the subpercent level.Unfortunately, for cool giants, each change made the interpolation error worse.
Interpolating model atmospheres boils down to independently interpolating four (for cool dwarfs) or five depth-dependent quantities over five stellar parameters.Which of these quantities and parameters are driving errors in cool stars?By replacing each depth-dependent quantity in a doubly interpolated atmosphere with the exact values, one at a time, we determined that the temperature at each layer is the largest driver of inaccuracy.Figure 5 demonstrates the effect of each stellar parameter on the interpolated temperatures at several layers.By varying each stellar parameter in turn, and directly examining the interpolated temperatures, we can identify the parameters that would benefit from finer sampling.While the temperatures at various optical depths are nearly linear as a function of log g, and [M/H], they vary sharply as a function of T eff , [α/M] and [C/M].This is because of the sensitivity of molecular formation to these parameters.For the case of [C/M], we also plot the number density of CO at the top of the atmosphere.It's clear that the sharp temperature changes at [C/M] ≈ 0.25 is driven by opacity from CO.We believe that other sharp transitions are similarly linked to formation of, and opacity from, specific molecules.
Finally, we note that the overwhelming challenges of model atmosphere interpolation for cool stars via standard techniques likely affects many analyses in the literature, as we are unaware of any past work documenting it.As shown in Figure 5, finer sampling in most parameters is required to enable accurate interpolation in this regime.In addition to denser sampling, additional abundance ratios will likely need to be included in model atmosphere grids to allow for accurate reconstruction of cool atmospheres.As indicated in Gustafsson et al. (2008), nitrogen abundances should be included as a free parameter.That work also indicates that the microturbulent velocity, which is especially impactful for cool, line blanketed stars, should also be allowed to vary.Adhoc calculation of model atmospheres may ultimately prove to be the only solution to fully address these problems.The latest versions of Korg include substantial upgrades to the chemical equilibrium solver.Korg now self-consistently determines the electron number density, n e , when solving for chemical equilibrium, rather than assuming the value provided by the model atmosphere.A configurable warning is raised when Korg produces a value substantially different from that in the model at-mosphere.The impact on the spectrum is nearly always undetectable, but allowing for n e to vary from the precise value in the model atmosphere allows Korg's chemical equilibrium solver to converge for almost all stellar parameters.Direct comparisons to the MARCS chemical equilibrium tables4 indicate that agreement is excellent except at the lowest temperatures, where small (≲ 10 −3 ) differences in Na and K ionization fraction drive moderate relative differences in the (very small) electron pressure.
In addition, Korg now supports polyatomic molecules, and includes partition functions for (the most abundant isotopologue of) each polyatomic molocule included in ExoMol5 version 20220926, except cis-P 2 H 2 and trans-P 2 H 2 (28 molecules total, listed in Table 1).In order to compute equilibrium constants, atomization (dissociation) energies are required for each molecule.These were calculated from the enthalpies of formation at 0 K of the molecule and its constituent atoms using experimental data from release 22 of the Computational Chemistry Comparison and Benchmark DataBase6 .Conventionally, nuclear spin degeneracy is not included in molecular energy level degeneracies in stellar astrophysics, and care should be taken that all data, including linelist log gf values and partition functions use the same convention.We obtained nuclear spin data from the Brookhaven National Laboratory "Wallet Card" ser- vice 7 and used it to convert the ExoMol partition functions from the "physics" convention (including the nuclear spin degeneracy) to the "astrophysics" convention.Finally, Korg now includes support for positively charged molecules.By default, those with parition function in Barklem & Collet (2016): H 2 +, He 2 +, C 2 +, N 2 +, O 2 +, Ne 2 +, P 2 +, S 2 +, HeH+, BeH+, CH+, NH+, OH+, HF+, NeH+, MgH+, AlH+, SiH+, PH+, SH+, HCl+, ZnH+, HBr+, CdH+, HgH+, CN+, CO+, NO+, NS+, BO+, SiO+, PO+, SO+, AsO+, and TaO+.

H I bf absorption and plasma effects
When atoms are embedded in a sufficiently hot and dense environment, their internal structure can no longer be dealt with in isolation.The Mihalas-Hummer-Däppen (MHD) formalism (Hummer & Miha-7 https://www.nndc.bnl.gov/nudat3/indxsigma.jsplas 1988) provides a probability that each level is dissolved into the continuum, and thus a correction to the partition function.The formalism was initially developed to provide an equation of state for stellar interiors, but by changing level populations (of, e.g.hydrogen), it also predicts changes to the monochromatic opacity.Plasma effects are typically subtle in stellar atmospheres, so the effect of MHD on the equation of state is negligible.Figure 6 shows the fractional correction to the neutral atomic hydrogen (the species most strongly affected by plasma effects) partition function evaluated at τ ros = 2.5, deep enough to have minimal impact on the spectrum.The largest corrections are for the coolest main-sequence stars, but all corrections are by factors of less than 4 × 10 −3 .Higher in the atmosphere, at depths that more strongly impact the emergent spectrum, the correction is smaller: 2 × 10 −4 at τ ros = 1 and 3 × 10 −5 at τ ros = 10 −2 .The corrections to level populations of hydrogen, however, are not negligible (see Figure 2  for both bound-bound (line) and bound-free (photoionization) opacity.Korg now includes the effects of the MHD formalism for hydrogen bound-free opacity, which has a particularly strong impact on the shape of the Balmer break.Due to level dissolution, photons that would ordinarily excite the atom to an upper energy level (i.e. they have energy less than the classical ionization energy -∼10.2eV from the n = 2 level) have a finite probability of ionizing the atom instead.This results in a "roundedoff" Balmer break.This calculation requires knowledge of the ionization cross section for energies lower than the classical ionization energy.We extrapolate the cross sections assuming that they are proportional to E −3 , where E is the photon energy.Using linear extrapolation instead results in differences in the cross section of up to a few percent.Applying the MHD formalism to the Lyman-α transition results in additional unphysical absorption far into the visible for any reasonable method  of extrapolating the photoionization cross section.Because of this, we do not apply the MHD formalism to the Lyman series.Wheeler et al. (2023) highlighted the fact that the application of MHD may over-attenuate hydrogen lines (see Section 3.3 for a related discussion).It has been noted previously that MHD dissolves outer, sparsely populated orbitals more eagerly that the OPAL (Rogers et al. 1996) equation of state (e.g.Iglesias & Rogers 1995).Nayfonov et al. (1999) presents a modified version of the MHD formalism that may address this in the stellar atmosphere regime, but further investigation is required.While we include the effects of MHD on hydrogen bound-bound opacity by default, we provide a switch to turn it off if desired.

Brackett lines: plasma effects and broadening mechanisms in the infrared
Spectral synthesis codes must treat hydrogen lines specially (if they are to model them accurately), because they are subject to resonant and linear Stark broadening.Support for Brackett series lines has been added to Korg, extending its applicability further into the infrared (out to about 2.3 µm, the Pfund series limit).Korg's implementation is adapted from the HLINOP hydrogen opacity code (Barklem & Piskunov 2003, 2016) and translated into Julia to facilitate automatic differentiation.HLINOP handles Brackett lines using the theory developed in Griem (1960), which extended the classic Holtsmark theory (Holtsmark 1919) to include the effects of ions, and Griem (1967), which added a correction for distant electrons.Comparisons between Synspec (Hubeny & Lanz 2011, 2017;Hubeny et al. 2021) and Turbospectrum (Plez et al. 1993;Plez 2012;Gerber et al. 2022) (used to generate model spectra for APOGEE in DR 17, Abdurro'uf et al. 2022) revealed that agreement on the Brackett lines is particularly poor, in contrast with other lines in the infrared, where synthesis codes tend to agree well (Wheeler et al. 2023).Since Korg's implementation was heavily inspired by HLINOP, the initial implementation produced Brackett lines which matched Turbospectrum's almost exactly, but the lines produced by Synspec differed significantly.Upon investigation, two independent factors were determined to be driving the disagreement between models: treatment of level dissolution, and oversights regarding line broadening mechanisms in the infrared.The first factor is simple: Turbospectrum includes the effects of the MHD formalism on hydrogen lines, while Synspec seems not to.The second factor requires slightly more explanation.
The top panel of Figure 7 shows the depth of formation for two hydrogen lines: H-α (visible), and Br-η (infrared).The core of H-α forms high in the atmosphere, at a Rosseland mean optical depth of 10 −3.5 or so.Because the monochomatic opacity is lower in the infrared than in the visible, the core of hydrogen lines in the infrared form much deeper, typically below the Roseland mean photosphere.This contrast has important implications for line broadening mechanisms-hydrogen lines that form deep in the atmosphere don't have Dopplerdominated cores.The second panel of Figure 7 shows the characteristic widths of hydrogen line profiles from the Doppler broadening, linear Stark broadening, and from other mechanisms (fundamental line widths and van der Waals broadening).At the depth of formation of H-α (marked with a dashed line), the Doppler broadening profile is the widest, so it sets the shape of the absorption profile resulting from convolving all three profiles.But at the depth of formation of Br-η, Stark broadening sets the shape of the line core.While none of the above constitutes novel theoretical insight, synthesis codes developed for application to the visible and near ultraviolet have generally overlooked the fact that hydrogen lines in the infrared need separate consideration.
The reason for the oversight is that numerically convolving broadening profiles is computationally expensive (in decades past, prohibitively so).As a consequence, all spectral synthesis codes commonly used for fitting data employ a variant of the same approximation: the line wings are obtained by adding all profiles, and the The characteristic profile widths of line broadening mechanisms as a function of depth.The dashed lines in each panel mark the Rosseland mean optical depth at which the core of each line forms.For line cores in the visible, which form high in the atmosphere, Doppler broadening dominates, but for line cores in the infrared, which form deep in the atmosphere, Stark broadening dominates.These calculations were done using a model atmosphere with T eff = 6000 K, log g = 3.5.
line core is set to the broadest one.This approximation works remarkably well, though problems arise at the boundary between the core and the wings, and when two of the broadening profiles have similar widths (this case arises for Brackett lines in cool stars, where the Stark and Doppler widths at the depth of formation can be comparable).For some codes, the issue is that the core is hard coded to be Doppler-dominated.For codes whose hydrogen lines descend from the Kurucz routines, the issue is more subtle.The Griem (1960Griem ( , 1967) ) theory employed provides separate profiles for Stark broadening in the quasistatic and impact limits.Profiles for each must be calculated using the number of perturbers (approximately) in each regime, then the profiles are con-volved to give the total Stark profile.These codes add the two profiles together, which gives a good approximation to the profiles in the wings, but results in a core with far too much opacity, and wavelength-integrated cross sections which are wrong by a factor of ∼ 2 when Stark broadening is dominant.Thankfully, this problem is easily remedied by using true numerical convolution of the broadening profiles.Figure 8 illustrates the effects of level dissolution and of broadening treatment on the Brackett series.It shows Br-η (the strongest hydrogen line in the APOGEE wavelength range), and Br-λ (a weaker line) for a few stellar parameters, as predicted by Synspec, Turbospectrum, and Korg.For each star, the probability that the upper level of each transition is not dissolved into the continuum, w, is plotted as a function of optical depth, with a gray band marking the depth of formation for the Brackett series.In some cases, MHD results in highly attenuated lines.Contrast Synspec and Korg without dissolution to Turbospectrum and Korg with dissolution in Br-λ in row 2 and in Br-η in row three.When level dissolution is weak (both lines in row 1, Br-η in row 2), the differences in predicted line profiles are due to differences in the broadening of the line core.
Figure 3.3 compares our improved Brackett profiles to a representative selection of Gaia FGK benchmark stars (Blanco-Cuaresma et al. 2014) observed by APOGEE.We show the Br-η lines computed using the naive and corrected profiles, and with level dissolution switched on and off.The lines are synthesized using the parameters from Blanco-Cuaresma et al. (2014), including nonspectroscopic T eff and log g.These comparisons clearly indicate that the corrected profiles significantly improve the accuracy of the synthesized spectra.Unfortunately, none of the benchmark stars provide a clear constraint on the effects of level dissolution, since including it doesn't consistently improve or degrade the fit.Comparisons for all Gaia FGK benchmark stars observed for APOGEE DR 17 are in Appendix A.
Finally, we note that ideally stellar spectral synthesis codes would treat the Brackett series using the model microfield method (Brissaud & Frisch 1971), as is often done for other hydrogen series via the tables of Stehlé & Hutcheon (1999).In lieu of this, Vidal-Cooper-Smith (VCS) theory (Smith et al. 1969;Voslamber 1969;Vidal et al. 1970Vidal et al. , 1971) ) as calculated by Lemke (1997) may be a promising alternative.

Linelist parsing
In addition to Vienna Atomic Line Database8 (VALD) linelists, Korg now supports a MOOG linelists with isotope information, "Kurucz"-format molecular linelists, and Turbospectrum linelists.When parsing linelists with isostopic information, Korg adjusts the log gf value of each line to account for isotopic abundances, which can be specified by the user (except in the case of preadjusted VALD linelists).

CONCLUSIONS
We have presented several updates to the Korg code for stellar spectral synthesis.The version of the code used in this work is 0.27.1.The new capabilities include routines for fitting observed spectra via either synthesis or equivalent widths.The equivalent width fitting routine was tested against Meléndez et al. (2014) and found to produce iron abundances with a significantly lower scatter than the original analysis.
We have described the model atmosphere interpolation method employed by Korg and shown that it introduces minimal error to synthesize spectra for stars with T eff > 4000 K. Stars with T eff < 4000 K have atmospheres that contain molecules in significant quantities, which greatly complicate the dependence of atmospheric thermodynamic quantities on stellar parameters.Flux errors introduced by interpolation are smaller in the infrared than the visible, because lines in the infrared form at deeper atmospheric layers where thermodynamic quantities are better interpolated.While it it possible to obtain sub-percent atmosphere interpolation error for cool dwarfs, finer sampling in stellar parameters is required if we are to make interpolation error truly negligible.Additionally, more parameters and abundances than are typically taken into account (e.g.microturbulence) are likely to impact molecular densities and opacities.It remains to be seen whether sufficiently-dense grids of cool star model atmospheres are feasible to generate and use; ad-hoc model atmosphere calculations may be necessary.
Korg's chemical equilibrium solver has been significantly updated to facilitate synthesis for cool stars.Support for charged and polyatomic molecules has been added, with state-of-the-art partition functions from the ExoMol project.The electron number density is now calculated when solving the chemical equilibrium equations, meaning that syntheses with abundances varying from the model atmosphere are more internally consistent.The vertical gray band marks the approximate depth of formation of the Brackett lines.These spectra were synthesized using the APOGEE DR 17 linelist and convolved to the instrument resolution (R = 22, 500).
We now apply the Mihalas-Hummer-Däppen formalism to neutral atomic hydrogen, allowing the outer orbitals of hydrogen to dissolve into the continuum in dense regimes.While this treatment (or something similar) is both physically motivated and necessary to reproduce a "rounded-off" Balmer break, there are shortcomings to the theory applied in the stellar atmosphere regime.
Finally, we highlight the fact that synthesis codes produce Brackett series lines that are in strong disagreement with each other.We show that the causes of the disparity are the treatment of level dissolution and that the fact that Stark broadening dominates hydrogen line cores in the infrared has been frequently overlooked.Comparisons to spectra of stars with non-spectroscopic values of T eff and log g indicate that our corrected Brackett profiles are in much better agreement with the data than those of other codes.investigation into H I bound-free absorption and plasma effects, and Andrew Saydjari (Harvard) for extensive stress testing of the code.

Figure 3 .
Figure 3.Estimated error introduced by interpolation as a function of T eff and log g for solar abundances.The color of each point shows the largest flux deviation at any wavelength for spectra synthesized with the APOGEE DR 17 linelist (Abdurro'uf et al. 2022) and convolved to R = 22, 500 (the approximate APOGEE resolving power).We note, however, that this plot is not sensitive to wavelength range, linelist, or R. For reference, the distribution of stars in APOGEE DR17 is plotted in the background.The colorbar has been truncated for clarity; errors are smaller than 10 − 6 for some stellar parameters.

Figure 4 .
Figure 4.Estimated error introduced by interpolation using the method tailored for cool dwarfs, described in text.Changes which improved interpolation for cool dwarfs make it worse for cool giants.(Compare the main sequence and red-giant branch to those in Figure 3.)

Figure 5 .
Figure5.The temperature at selected upper layers in interpolated (using the resampled cubic method) model atmospheres (10 −4 ≲ τ5000 ≲ 10 −1 ) in exact (black) and interpolated (dashed red) model atmospheres.The interpolated atmospheres are always interpolated in all five stellar parameters.Each panel shows how the values change as one of the parameters varied and the others stay fixed.The fixed values (T eff = 3800 K, log g = 0.0, solar abundances), are marked with vertical lines.In the last panel, the blue line indicates the number density of CO at the top of the model atmosphere, computed using the uninterpolated atmospheres only.Sharp changes in temperature as a function of stellar parameters driven by molecular formation make model atmosphere interpolation very challenging for cool giants.

Figure 6 .
Figure 6.Fractional correction to the hydrogen partition function, UH , deep in the stellar atmosphere according to the MHD formalism across the Kiel diagram.Stellar atmospheres are produced using Korg's interpolation of the SDSS model atmosphere grid for solar-composition stars.Corrections are less than 4 × 10 −3 for all stars and much smaller than that at higher layers.

Figure 7 .
Figure 7. Top: the Rosseland mean optical depth of formation of H-α and Br-η as a function of wavelength.Lines in the infrared form much deeper in the stellar atmosphere, where different broadening mechanisms dominate.Bottom:The characteristic profile widths of line broadening mechanisms as a function of depth.The dashed lines in each panel mark the Rosseland mean optical depth at which the core of each line forms.For line cores in the visible, which form high in the atmosphere, Doppler broadening dominates, but for line cores in the infrared, which form deep in the atmosphere, Stark broadening dominates.These calculations were done using a model atmosphere with T eff = 6000 K, log g = 3.5.

Figure 8 .
Figure8.Predicted lines of Br-η and Br-λ (n = 4 to n = 11 and n = 15, respectively), two lines in the APOGEE wavelength range.Agreement between synthesis codes is poor because models don't properly account for line formation in the infrared and because of differing treatments of level dissolution.Each row represents a different set of stellar parameters (T eff and log g are indicated in the first column, and all models are computed with [M/H] = −1.5).The third column shows the probability that the n = 11 and n = 15 orbitals survive (are not dissolved into the continuum) as a function of optical depth at 5000 Å, τ5000.The vertical gray band marks the approximate depth of formation of the Brackett lines.These spectra were synthesized using the APOGEE DR 17 linelist and convolved to the instrument resolution (R = 22, 500).

Figure 9 .
Figure9.Comparisons of synthesized Br-η profiles for four representative Gaia FGK benchmark stars.The fit is much improved using the corrected profiles, though the effect of the MHD dissolution formalism is to subtle to constrain using these data.
Meléndez et al. (2014)Fe/H] abundances of 18 Sco (relative to the Sun) computed by Korg with equivalent widths and stellar parameters fromMeléndez et al. (2014).The mean abundance is the same found byMeléndez et al. (2014), and the scatter in line abundances from Korg is nearly half that ofMeléndez et al. (2014)(0.0056dex compared to 0.010 dex).

Table 1 .
in Wheeler et al. 2023), which has implications Polyatomic molecules from the ExoMol project which have been added to Korg