© ESO 2020

1 Introduction

According to Standard Big Bang Nucleosynthesis (SBBN), ⁷Li together with hydrogen and helium, is one of the few elements created in the primordial fireball. To trace the primordial chemical composition of the Universe, one of the best diagnostics is represented by old Population II stars, characterized by effective temperatures T_eff ≳ 5900 K. Such objects are metal-poor and their convection zones are shallow enough to preserve their initial lithium content. More than 30 yr ago, Spite & Spite (1982) analyzed a sample of unevolved metal-poor (MP) stars (− 2.8 ≤[Fe∕H] ≤−2.0) with temperatures in the range 5700 ≤ T_eff ≤ 6300 K and discovered that irrespective of effective temperature, gravity, and metallicity, all analyzed stars had essentially the same Li abundance of A(Li)¹ ≈ 2.2 that defines the Spite plateau. The invariable behavior of lithium at low [Fe∕H] is different from what occurs for other elements, whose abundances generally drop with declining metallicity. The Spite plateau has been interpreted as the primordial lithium abundance, produced during the SBBN, and preserved in the atmosphere of these stars due to their shallow convection zones. The fact that this constant abundance value is three times smaller than the primordial Li abundance inferred from the measurements of the cosmic microwave background (CMB) by WMAP (Spergel et al. 2007) and Planck (Planck Collaboration XVI 2014) poses a severe problem for the interpretation of the Spite plateau as representative of the primordial Li abundance. This puzzling discrepancy challenges our present understanding of Galactic chemical evolution and stellar Li depletion.

Recent studies (e.g., Korn et al. 2006, 2007; Nordlander et al. 2012) explain this discrepancy by invoking a combination of diffusion processes and some form of turbulent mixing operating at the bottom of the outer convective zone in turn-off stars belonging to the Spite plateau, e.g. in globular cluster NCG 6397. As discussed in Richard et al. (2005), this extra ingredient is needed to meet the constraints imposed by the observations of ⁷Li in old halo stars: a thin and flat plateau can only be retained in the presence of diffusion if extra mixing is introduced into the models. This process rapidly transports the surface Li down to the hotter inner region, where it is readily destroyed, resulting in a Li-depleted convection zone. Since there is no parameter-free physical description of the turbulent mixing, it is necessary to introduceit into the stellar evolution models in an ad-hoc manner, defining a turbulent diffusion coefficient D_T, which is a parametric function of density and the atomic diffusion coefficient of helium at a reference temperature T₀ (; see Richard et al. 2005). High-quality A(Li) measurements are fundamental for understanding whether or not atomic diffusion is at work in metal-poor stellar atmospheresand if so, to what extent. In their Fig. 5, Nordlander et al. (2012) show the behavior of different diffusion models in comparison with observed stars in globular cluster NCG 6397. Their measured 1D non-local thermodynamical equilibrium (NLTE) lithium abundances for the two metal-poor stars belonging to the Spite plateau (at T_eff ≈ 5900 and 6300 K) can be fitted, within error bars, with one of their fine-tuned diffusion models. The predicted primordial lithium abundance A(Li) = 2.57 ± 0.1 is compatible with the primordial Li abundance inferred from Planck data. The same model, however, fails to reproduce the red-giant branch (RGB) stars at lower temperatures, probably due to a source of extra mixing that was not considered in these calculations (Charbonnel 1995).

The same globular cluster was also analyzed by González Hernández et al. (2009), who measured the 3D NLTE lithium abundance in a sample of subgiants (SG) and main-sequence (MS) stars, aimed at investigating the cosmological Li discrepancy in the framework of diffusion models by tracking the lithium abundance trend with evolutionary phase. As shown in their Fig. 2, they found SG stars to have a higher lithium abundance than MS objects and in both cases, a decreasing trend of A(Li) with T_eff with similar slopes was observed. Although no conclusive explanation of such trends was provided, they compared the results with different diffusion models but were not able to fit satisfactorily the A(Li) versus T_eff trend of their targets, starting from an initial (primordial) Li abundance in agreement with the Planck observations. The only model able to provide a decreasing T_eff trend, with abundances of SG stars larger than MS stars, failed quantitatively because the predicted A(Li) for their warmest stars was about 0.05 dex lower than observed, such that the overall slope of A(Li) with T_eff was not properly reproduced. Similar difficulties were encountered in the work by Gruyters et al. (2016), again in the cluster NGC 6397 while fitting their abundances of lithium and other elements (Fe, Ca) as a function of T_eff with diffusion models. The authors advise caution about the size of the abundance trends predicted by the models, which is set by the relative abundance difference between the hottest and the coolest stars of the target sample.

These examples demonstrate how crucial it is to derive lithium abundances with a precision of better than 0.1 dex, needed to distinguish between different diffusion models and to better understand the underlying physics. Furthermore, this would help in constraining the efficiency of atomic diffusion and additional turbulent mixing processes, which at the current stage represents a plausible scenario for explaining the gap in lithium abundance as measured in unevolved MP stars and inferred from the CMB measurements by WMAP and Planck.

Another open issue is related to the question of whether the Spite plateau presents some slope with the metallicity or if it is actually a flat distribution of A(Li) (see e.g., Thorburn 1994; Norris et al. 1994; Ryan et al. 1996, 1999; Asplund et al. 2006; Meléndez et al. 2010). The plateau, in fact, presents a flat trend with some dispersion in lithium abundance between [Fe∕H] = −2.2 and − 2.5, as originally found by Spite & Spite (1982), and the observed scatter can be as low as 0.033 dex, as found in Ryan et al. (1999). At metallicities lower than [Fe∕H] ≤−2.5, instead, a larger dispersion in the Li abundances is observed with a trend of decreasing A(Li) with decreasing [Fe∕H] (Bonifacio et al. 2007; Aoki et al. 2009; Sbordone et al. 2010). At present, it is not clear from neither an observational nor a theoretical perspective, why most stars in this metallicity range show lithium abundances below the Spite plateau. More recently, Bonifacio et al. (2018) revisited this so-called “meltdown” effect of the Spite plateau, finding some stars at [Fe∕H] < −3.0 with A(Li) below the constant value around ≈ 2.2 dex (as derived by Sbordone et al. 2010)but also two stars at [Fe∕H] ≈−4.0 with a lithium content that is fully compatible with the usual flat distribution of lithium in MP stars (the two star symbols in Fig. 2 of Bonifacio et al. 2018 at A(Li) = 2.2 and [Fe∕H] ≈−4.0). The presence of these two stars is exceptional since, as pointed out by Matsuno et al. (2017), no star in the literature was known to have a Li abundance comparable to the Spite plateau below [Fe∕H] of − 3.5, except for the primary of the double-lined binary system CS 22 876-032 (González Hernández et al. 2008). Very recently, an even more extreme case was found: Aguado et al. (2019) discovered that SDSS J0023+0307, one of the most iron-poor stars known ([Fe/H] ≲ 6.1), lies on the Spite plateau with A(Li) = 2.02 ± 0.05.

In order to provide a better observational basis for meaningful comparisons with stellar depletion models and to shed more light on the open questions regarding the Spite plateau, precise Li abundances for a large sample of MP stars are needed, accounting for both 3D and NLTE effects that may affect significantly the strength of the lithium line in metal-poor stars and, hence, their estimated lithium content. Highly accurate lithium abundances are also required in statistical studies trying to detect significant differences in the lithium content of the atmospheres of planet host stars and stars without known planets (e.g., Ghezzi et al. 2009; Israelian et al. 2009; Baumann et al. 2010; Pavlenko et al. 2018).

In this work, we present extended fitting functions to derive the Li abundance of solar-type stars of different metallicity, taking into account both hydrodynamical effects (3D) and departures from thermodynamical equilibrium (NLTE). These fitting functions allow to derive (i) the 3D and 1D NLTE abundance correction; (ii) the 3D and 1D NLTE equivalent width corresponding to a given Li abundance A(Li); (iii) the NLTE lithium abundance A(Li)_3DNLTE and A(Li)_1DNLTE from measured equivalent widths. Our approach is similar to that of Harutyunyan et al. (2018) who additionally provide 3D NLTE corrections for the lithium isotopic ratio, albeit for a smaller range of stellar parameters. The present work extends and supersedes the fitting functions presented in Sbordone et al. (2010).

2 Stellar atmospheres and spectrum synthesis

2.1 3D CO5BOLD and 1D LHD model atmospheres

For this work, we adopted a sub-set of the CIFIST grid of 3D hydro-dynamical CO5BOLD model atmospheres (Ludwig et al. 2009)together with their corresponding 1D LHD models. These 1D models employ the same stellar parameters, micro-physics, chemical composition, and radiative transfer scheme as a CO5BOLD model, facilitating an unbiased 3D–1D differential comparison (Caffau et al. 2007). Our model grid is shown in Fig. 1 in the log g− T_eff plane that can be seen as a pseudo Hertzsprung-Russell (HR) diagram. The grid includes 70 3D models focusing on main-sequence and turn-off coolstars, represented by the following stellar parameters:

• T_eff: 5000, 5500, 5900, 6300 and 6500 K,

• log g: 3.5, 4.0 and 4.5 [cm s⁻²],

• [Fe∕H]: 0.00, −0.5, −1.0, −2.0 and −3.0 dex.

The region around T_eff = 6500 K and log g=3.5 (blue filled squares in Fig. 1) is only scarcely populated in the (observed) HR diagram, and for this reason models with such parameters are not available in the CIFIST grid. However, since a rectangular T_eff – log g grid is more practical for defining a homogeneous fitting function, we decided to provide approximate results also for this missing combination of parameters by extrapolation of the lithium abundance corrections derived from the available atmospheres.

In Table A.1, the main properties of the grid of selected models are listed. The effective temperature reported in the fourth column of Table A.1 is the average T_eff over the sub-sample of N_snap selected snapshots (where N_snap is given in the third column) for each CO5BOLD simulation. Note that these T_eff deviate slightly from the nominal effective temperatures T_nom defining the rectangular grid depicted in Fig. 1 (where T_nom=5000, 5500, 5900, 6300 and 6500 K for all log g and [Fe∕H]). The fitting function will be developed on this fixed temperature grid after interpolation of the relevant quantities from the true T_eff to T_nom.

For [Fe∕H] = 0.0 and − 2.0, we have in addition computed independent, refined grids of 1D LHD models, covering the same range of gravities (log g = 3.5, 4.0, 4.5) and effective temperatures (5000 ≤ T_eff ≤ 6500 K), but with a smaller and equidistant step of ΔT_eff = 100 K. These auxiliary grids of 1D LHD models are useful for estimating the interpolation errors of the fitting functions developed in Sect. 4.

Fig. 1 Overview of the grid of CO5BOLD 3D model atmospheres adopted for the computation of the fitting functions. The green filled squares mark the models that have been used in this work (five different metallicities, see legend at upper left), the grey squares indicate existing models that are part of the CIFIST grid but have not been considered in the present context. The blue squares highlight the missing models for which abundance corrections were derived by extrapolation (see Sect. 3.2). The solar 3D CO5BOLD model is indicated by the Sun symbol.

2.2 Library of 3D and 1D NLTE synthetic spectra

To ensure that our fitting functions cover a sufficiently wide range of lithium abundances, we consider eleven values of A(Li), from 1.0 to 3.0 with a step of 0.2 dex. For each abundance, we calculate individual NLTE departure coefficients for all snapshots of each 3D model and for their associated one-dimensional LHD counterparts, using our NLTE3D code (see, e.g., Steffen et al. 2015; Mott et al. 2017). We adopted the same lithium model atom as described in Mott et al. (2017), which is more complete than the one used by Sbordone et al. (2010), who, moreover, calculated the NLTE departure coefficients only for one representative value of A(Li) = 2.2 and assumed them to be valid for all considered lithium abundances in the range 0.9 ≤ A(Li) ≤ 3.3.

Our approach of calculating a separate set of departure coefficients for each abundance constitutes a more reliable treatment of the NLTE effects over the full range of line strengths than the simplified approach adopted by Sbordone et al. (2010). In this sense, their fitting functions are superseded by the updated and extended version provided in the present work.

Once the NLTE departure coefficients are available, they are fed into the NLTE spectral synthesis code Linfor3D² to compute a library of synthetic lithium line profiles. Our library consists of 770 3D NLTE lithium line profiles in the wavelength range λ = [670.69 − 670.87] nm (14 log g –T_eff combinations × 5 metallicities × 11 Li abundances) where each spectrum is the average over the selected snapshots. For each of these 3D NLTE spectra, we have computed corresponding 1D NLTE and 1D LTE spectra for comparison. In total, the library of synthetic spectra then comprises 3 × 770 = 2310 spectra. In addition, we have computed a LTE curve-of growth (CoG) for each set of stellar parameters (T_eff, log g, [Fe∕H]) which allows us to interpolate the 1D LTE equivalent width (EW_1DLTE) of the lithium line over a wide range of A(Li) between 0.0 and 4.5 dex. In exactly the same way, we calculated a total of 1056 additional 1D NLTE line profiles (48 log g–T_eff combinations × 2 metallicities × 11 Li abundances) plus the corresponding LTE CoGs from the auxiliary (refined) grids of 1D LHD models, available for [Fe∕H] = 0.0 and − 2.0 only.

2.2.1 Microturbulence in 1D LHD models

In the 1D models, the small-scale velocity fields in the stellar photosphere is represented by the so-called microturbulence parameter ξ_micro. As usual, ξ_micro describes a depth-independent isotropic Gaussian velocity distribution with a dispersion of v_rms = ξ_micro.

In the 3D models, on the other hand, there is no need to introduce such artificial fudge parameter since the full velocity field is obtained naturally from the hydro-dynamical treatment of the convective motions in each CO5BOLD model atmosphere.

For a differential 3D-1D comparison, the microturbulence parameter of the 1D LHD model associated with a given 3D model should be adjusted such that the resulting small-scale velocity field matches (as closely as possible) the properties of the 3D hydrodynamical velocity field. Such a calibration has been presented by Dutra-Ferreira et al. (2016) who provide the following (theoretical) relation for ξ_micro as a function of effective temperature and surface gravity: (1)

where X = T_eff − 5500 [K] and Y = log g − 4.0.

This formula is the result of an analysis of 30 Fe I plus 7 Fe II lines for which synthetic spectra have been computed adopting 3D CO5BOLD and 1D LHD model atmospheres of solar metallicity (we refer to this paper for more details). Equation (1) constitutes the most appropriate description of the ξ_micro parameter available in the context of 3D–1D comparisons. Even though formally valid for solar metallicity only, we use this relation to assign ξ_micro to all LHD models used in the present study. This may be justified by empirical evidence that ξ_micro depends only weakly on [Fe∕H] (Sitnova et al. 2015; Mashonkina et al. 2017). The last column of Table A.1 lists the values of ξ_micro computed with Eq. (1).

2.2.2 Line list

The lithium resonance doublet at λ670.8 nm is represented by 12 hyper-fine structure components, six for each isotope, extracted from the original list of Kurucz (1995), as explained in Mott et al. (2017). For the spectrum synthesis we assume a solar system isotopic ratio of ⁶Li∕⁷Li = 8.2% (e.g., Lodders 2003).

Although our grid includes models with solar metallicity for which contributions to the equivalent width coming from some blend lines of different elements may be present, we excluded such features in the spectral synthesis, maintaining our focus on targets where these other atomic and molecular species have a negligible effect on the equivalent width (EW) of the lithium resonance doublet. See, however, Sect. 6.3, where we discuss the applicability of our abundance corrections in the presence of blend lines.

Once the grids of 3D NLTE theoretical spectra and 1D (both NLTE and LTE) CoGs for the λ670.8 nm lithium line have been completed, we initiated the procedure for deriving the abundance corrections for each model as a function of A(Li) and their transformation into the functional fits. These are the result of several steps that are described in the next sections.

3 3D NLTE and 1D NLTE lithium abundance corrections

3.1 Definition

When dealing with the concept of abundance corrections, it is important to define properly how such corrections are derived since they can be calculated in multiple ways. Considering the 1D NLTE case, the corrections, Δ_1DNLTE-1DLTE, are usually defined as a function of the 1D LTE abundance, , such that: (2)

To derive this correction for a given (which is usually the result of a standard 1D LTE analysis), one first has to compute the 1D LTE equivalent width associated with this LTE lithium abundance, EW_1DLTE. The next step is the search for the 1D NLTE lithium abundance, that provides the same EW_1DLTE value. Note that this step requires the computation of a 1D NLTE curve-of-growth. The sought 1D NLTE abundance correction (Δ_1DNLTE-1DLTE) is then defined as the difference of the two abundances. These corrections can be applied to standard 1D LTE abundances to correct them for non-LTE effects.

A completely analogous procedure could be applied to derive the 3D NLTE corrections Δ_3DNLTE-1DLTE, again as a function of the 1D LTE abundance , allowing to correct the 1D LTE abundance for the combined 3D plus NLTE effect: (3)

The “inverse” abundance correction is defined as the correction that has to be added to a given 3D NLTE abundance in order to obtain the corresponding 1D LTE abundance. It has the opposite sign but in general (for partly saturated lines) a different magnitude than the commonly used “forward” correction Δ_3DNLTE-1DLTE. While it is technically more straightforward to derive the “inverse” corrections from the theoretical models, they are of less practical interest and we calculate them only as an intermediate step. Eventually, we provide the “forward” corrections Δ_3DNLTE-1DLTE as a function of the 1D LTE lithium abundance that is typically derived by analyzing an observed stellar spectrum.

3.2 Computation of the lithium abundance corrections

The abundance corrections defined in Sect. 3.1 are derived by the following procedure.

Given our fixed set of lithium abundances (1.0 ≤ A(Li) ≤ 3.0), we computethe equivalent width of the Li I doublet from the 3D models in full 3D NLTE (EW_3DNLTE), and the LTE and NLTE EWs from the associated 1D LHD models (EW_1DLTE and EW_1DNLTE, respectively) by numerical integration of the synthetic lithium line profiles over the wavelength range Δλ = [670.69−670.87] nm. These data define three different curves-of-growth for each set of stellar parameters, corresponding to 3D NLTE, 1D LTE, and 1D NLTE treatment, respectively.

As an illustration, we show in Fig. 2 the 3D NLTE CoG (black) and the 1D LTE CoG (blue) for the models with T_eff = 5160 K, log g = 4.0 and [Fe∕H] = − 3.0. Filled dots mark the computed EWs at the prescribed grid of lithium abundances, continuous lines show the CoGs obtained by piecewise cubic interpolation. According to the above definition, the 3D NLTE abundance correction Δ_3DNLTE-1DLTE corresponds to the horizontal distance between the two CoGs. For a given 1D LTE abundance (e.g., A(Li) = 2.8), Δ_3DNLTE-1DLTE is negative and its magnitude is indicated by the length of the blue horizontal arrow, which can be easily obtained by interpolation.

For completeness, Fig. 2 also shows the “inverse” abundance correction , indicated bythe black horizontal arrow. This is the correction that has to be applied to a given 3D NLTE abundance in order to obtain the corresponding 1D LTE abundance. As mentioned above, it is not considered here as it is not of practical interest.

We derived the 1D NLTE corrections Δ_1DNLTE-1DLTE following the same method, replacing the EW_3DNLTE by the 1D NLTE equivalent widths EW_1DNLTE, computed with the microturbulence parameter that is listed in Table A.1 for each model atmosphere. The obtained grids of 3D NLTE and 1D NLTE abundance corrections have four dimensions: T_eff, log g, [Fe∕H], A(Li). At this stage, however, the grids need some further manipulation, before we are able to represent them by analytical expressions.

The first issue is that the corrections are based on equivalent widths that are computed for the real effective temperatures that depart slightly from the fixed nominal temperature grid T_nom = 5000, 5500, 5900, 6300, and 6500 K. This means that the temperature grid is slightly different for each (log g, [Fe∕H]) combination, preventing a straightforward comparison of abundance corrections for models with different surface gravity or metallicity at fixed T_eff. This difficulty is eliminated by interpolating the corrections from the actual T_eff grid to the fixed T_nom grid, separately for each (log g, [Fe∕H]) combination and for each of the 11 values of A(Li). Here we employ again a piecewise cubic interpolation, in some cases involving some very slight (uncritical) extrapolation. After this step, the grids of 3D NLTE and 1D NLTE abundance corrections are represented on a four-dimensional rectangular grid.

The second issue concerns the missing corrections for surface gravity log g = 3.5 and T_eff = 6500 K. Due to the lack of the corresponding model atmospheres in the CIFIST grid (blue squares in Fig. 1), we decided to find these values by means of extrapolation from the available corrections. For given [Fe∕H] and A(Li), we compute the missing corrections as the average of the results obtained from two orthogonal linear extrapolations: (i) using the available corrections at fixed log g = 3.5 and T_nom = (5900, 6300) K, and (ii) using the corrections at fixed T_nom = 6500 K and log g = (4.0, 4.5).

Even if we could have generated the missing 1D LHD models and the corresponding 1D LTE and NLTE spectra, we decided to compute the missing 1D NLTE corrections by extrapolation in exactly the same way as in the 3D case described above. Comparison with the independent results from the auxiliary (refined) grid of 1D LHD models for [Fe∕H] = 0.0 and − 2.0 (see Sect. 2.1) reveals that the extrapolation works perfectly for log g = 4.5, but may introduce an error of the order of − 0.01 dex for log g = 3.5 and 4.0.

After these preparatory steps, our rectangular grids of Δ_3DNLTE-1DLTE and Δ_1DNLTE-1DLTE corrections are complete. An overview of the corrections as a function of the grid parameters A(Li), T_eff, log g, and [Fe∕H] is presented in Fig. 3. The dependence of the 3D and 1D corrections on the stellar parameters and the lithium abundance is qualitatively similar. For weak lines (low A(Li)), they are generally positive, while they tend to be negative for the highest Li abundance. The corrections are most sensitive to the line strength for the coolest models. They are most negative for the lowest metallicities, reaching values of at least − 0.3 dex at [Fe∕H] ≤ − 2.0, T_eff = 5000 K, A(Li) = 3.0.

Fig. 2 Two different definitions of abundance corrections for the example of a strong Li line on the non-linear part of the CoG. The correction is indicated by the black horizontal arrow starting on the 3D NLTE CoG (black) at A(Li) = 2.8 and ending at a higher abundance given by the intersection with the 1D LTE CoG (blue). Conversely, the correction Δ3DNLTE-1DLTE is given by the blue horizontal arrow starting on the 1D LTE CoG (blue) at the same abundance A(Li) = 2.8 and extending to the 3D NLTE CoG (black).

4 Analytical approximations

4.1 Lithium abundance corrections

In the following, we derive a first parametric equation (hereafter FFI_3DNLTE, FFI_1DNLTE) that approximates the Δ_3DNLTE-1DLTE and Δ_1DNLTE-1DLTE corrections, respectively, as a function of T_eff, log g, [Fe∕H], and A(Li).

An effective way to achieve this is to perform, as a first step, a 2D-fit of the corrections in the T_nom-A(Li) plane for each [Fe∕H] and log g, separately. For this purpose, we utilized the IDL procedure SFIT which allows to perform a surface fit to a set of 2-dimensional data (in our case the Δ values given on a rectangular grid as a function of T_nom and A(Li)). SFIT returns an array of coefficients, the number of which depends on the chosen degree of the polynomial function used in the fitting.

Obviously, the goal is to achieve a reasonable fit of the corrections with a relatively low number of coefficients. After several tests we found that our corrections were satisfactorily fitted by a polynomial function with a maximum degree equal to 4, expressed as: (4)

where y = A(Li)_1DLTE (1.0≤ y ≤ 3.0) and z = (T_nom∕1000 − 5) (0 ≤ z ≤ 1.5). Such normalized coordinates make the 15 coefficients c_jk dimensionless and of similar magnitude.

Lowering the degree of the polynomial function, we noticed a significant degradation of the fitting quality for some metallicities or for some surface gravities. As a consequence, we conclude that in order to fit homogeneously the whole set of corrections with a unique polynomial form, the optimal choice is a maximum degree equal to 4.

Following this strategy, we obtain a set of 15 numerical coefficients for each given metallicity [Fe∕H] and gravity log g that defines the fitting surface of the corrections in the T_nom-A(Li) plane.

Noticing that the gravity dependence of the corrections is rather weak, we represent the gravity dependence of the coefficients c_jk by a linear function in the normalized surface gravity x = log g − 3.5 (0 ≤ x ≤ 1) as: (5)

The coefficients c_ijk are obtained from the c_jk given at the three gravities (x = 0.0, 0.5, 1.0) by standard linear regression, separately for each metallicity [Fe∕H].

For each particular metallicity value of our grid, we thus obtain the following final fitting function that describes the 165 NLTE lithium abundance corrections as a function of T_eff, log g, and A(Li) through 30 numerical coefficients: (6)

The list of 30 coefficients c_ijk are given, for each metallicity covered by our grid ([Fe∕H] = 0.0, −0.5, −1.0, −2.0 and −3.0 dex), in Tables B.1 (Δ_3DNLTE-1DLTE) and B.2 (Δ_1DNLTE-1DLTE). We decided not to include the metallicity as a fourth parameter of our fitting function since the trend of the corrections as a function of [Fe∕H] is not as simple as for the log g dependence. Finding a global fitting function able to fit our data for each combination of T_eff, log g, [Fe∕H], and A(Li) simultaneously would have required the inclusion of an impractical number of coefficients.

When applying this fitting function to metallicities other than those that are tabulated, a further interpolation of the obtained Δ values to the desired [Fe∕H] is necessary. We do not expect that this final operation would compromise the precision of the derived abundance corrections.

To facilitate the practical application of the analytical fitting function, we provide a Python script that computes the corrections as a function of the four parameters T_eff, log g, [Fe∕H], and A(Li), using the coefficients listed in Tables B.1 and B.2, respectively, according to Eq. (6), and subsequently performing a piecewise cubic monotonic interpolation (see Steffen 1990) in [Fe∕H]. An overview of the lithium abundance corrections obtained with the fitting functions FFI_3DNLTE and FFI_1DNLTE is shown in Figs. 4 and 5, respectively, as a series of contour plots in the T_eff – A(Li) plane, separately for each combination of [Fe∕H] and log g.

To check how well our fitting functions are able to represent the Δ_3DNLTE-1DLTE and Δ_1DNLTE-1DLTE input data, we computed for each combination of [Fe∕H], log g, T_eff, and A(Li) the difference between the abundance correction predicted by the fitting function and the value of the original input data. Figure B.1 shows these differences at the 4D grid points, both for the 3D NLTE and the 1D NLTE corrections.

For solar metallicity, the residuals δ_FFI fall in the range |δ_FFI|≲0.01, with a rms mean of . For all subsolar metallicities, the residuals are slightly larger, lying in the range |δ_FFI|≲0.02, with a typical rms mean of 0.005. With very few exceptions, the maximum residuals occur at T_eff = 5000 K for the extreme lithium abundances, A(Li) = 1.0 or 3.0.

In the 1D case, we can additionally evaluate the interpolation errors of FFI_1DNLTE for [Fe∕H] = 0.0 and − 2.0 by comparison with the corrections derived independently from the auxiliary grid of 1D LHD models for a refined sample of effective temperatures. We find that the interpolation errors at the intermediate temperatures have about the same amplitude as the fitting function errors at the nodes of the original grid, with |δ_FFI|≲0.016 dex in all cases, as shown in full detail in Fig. B.2. At an arbitrary point in the multidimensional parameter space (T_eff, log g, A(Li)), the interpolation error may be somewhat larger, due to the additional (unknown) interpolation errors in log g and A(Li). However, we can assume that the latter two interpolation errors are small compared to the interpolation error in T_eff, due to the weak gravity dependence of the corrections and the fine spacing of the original A(Li) grid, respectively, such that the interpolation errors shown in Fig. B.2 can still be considered as a valid order-of-magnitude estimate of the total interpolation error.

Taken together, the error analysis indicates an overall good fit of tabulated abundance corrections by our fitting function FFI_1DNLTE, and by implication FFI_3DNLTE, across the whole T_eff –log g–A(Li) space. As an additional benefit, the fitting functions smooth out small numerical artifacts in the raw data.

Fig. 3 Overview of the Li abundance corrections Δ3DNLTE-1DLTE (top) and Δ1DNLTE-1DLTE (bottom) as evaluated on the complete 4D parameter grid. For each metallicity, a total of 165 points indicates the corrections as a function of an arbitrary spectrum number that runs over A(Li), Teff, and log g. Large (left group), intermediate (middle group), and small dots (right group) correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right (5000 K: red, 5500 K: black, 5900 K: gray, 6300 K: black, 6500 K: blue). Each effective temperature is represented by 11 dots indicating the results for 1.0 ≤ A(Li) ≤ 3.0.

Fig. 4 Contour plots of the 3D NLTE abundance corrections in the Teff–A(Li)1DLTE plane. Each panel shows the corrections computed with FFI3DNLTE (Eq. (6)) for all combinations of surface gravity and metallicity as given in the title of each plot.

Fig. 5 Contour plots of the 1D NLTE abundance corrections in the Teff–A(Li)1DLTE plane. Each panel shows the corrections computed with FFI1DNLTE (Eq. (6)) for all combinations of surface gravity and metallicity as given in the title of each plot.

4.2 From A(Li) to EW

In this section we describe the derivation of a second parametric equation (hereafter FFII_3DNLTE, FFII_1DNLTE) which allows the computation of the 3D NLTE and 1D NLTE equivalent width values for the input parameters T_eff, log g, [Fe∕H], and lithium abundance A(Li).

The principal ingredient is obviously represented by the equivalent width values of the lithium doublet for each one of the 70 model atmospheres in our grid as computed for each one of the 11 lithium abundances in the range A(Li) = [1.0 .. 3.0]. For the 1D case only, the microturbulence parameter ξ_micro of each model has been set to the value calculated by means of Eq. (1). Otherwise, the derivation of the functional fit described in the following is the same for the 3D NLTE and the 1D NLTE case.

In a first step, the EWs obtained from the model atmospheres (at given log g, [Fe∕H], and A(Li)) were interpolated from the actual T_eff values to the fixed T_nom grid. In a second step, the resulting rectangular grid of EW values is completed by estimating the missing EW values at T_nom=6500 K and log g=3.5 for each metallicity by means of linear extrapolation as described in Sect. 3.2.

For each log g – [Fe∕H] combination, we thus obtain a complete, rectangular grid of EW values as a function of temperature and lithium abundance. The derivation of the fitting function FFII is completely analogous to the procedure to derive the fitting function FFI for the abundance corrections as described in Sect. 4.1. First we use the IDL procedure SFIT to compute the coefficients of the two-dimensional polynomial function of maximum degree four providing the best surface fit to the (logarithmic) EWs in the y – z plane, separately for each of the 15 log g – [Fe∕H] sets. Then the log g-dependence of the coefficients is represented by a linear function in x = log g − 3.5. In this way, we obtain for each of the five metallicities a set of 30 dimensionless coefficients w_ijk that define the fitting function FFII for the EWs. It has the same form as Eq. (6): (7)

where EW is the equivalent width in mÅ, and as before x = log g − 3.5, y = A(Li), z = (T_nom∕1000 − 5). The coefficients w_ijk are given in Tables B.3 and B.4 for the 3D NLTE and 1D NLTE case, respectively.

The comparison plot shown in Fig. 6 demonstrates how well our function fits the 3D NLTE input data. Clearly, the FFII fitting polynomial of maximum degree four (continuous red lines) ensures a very good agreement between the input EW (filled dots, color coded according to T_nom) across the full parameter space covered by our grid. The fit of the 1D NLTE EWs is even slightly better. For both 3D and 1D, the relative fitting function error is |ΔEW/EW| ≲ 2% with few exceptions, while the rms mean fitting error on the 4D parameter grid is less than 1%. The maximum global fitting function error is ΔEW/EW = + 3.9% for FFII_3DNLTE ([Fe∕H] = − 3.0, log g = 4.0, T_eff = 5000 K, A(Li) = 1.0) and − 3.4 % for FFII_1DNLTE ([Fe∕H] = − 1.0, log g = 4.0, T_eff = 5900 K, A(Li) = 1.0). Interpolation errors are of similar size. For further details see Figs. B.3 and B.4.

For practical application of FFII, we provide a Python script that computes the EWs as a function of the four parameters T_eff, log g, [Fe∕H], and A(Li), using the coefficients listed in Tables B.3 and B.4, respectively, according to Eq. (7), followed by a piecewise cubic monotonic interpolation in [Fe∕H].

Fig. 6 Comparison between the input data (filled dots, color-coded according to effective temperatures as shown in the legend of the upper left panel) and the resulting fitting polynomial FFII3DNLTE (Eq. (7), red continuous line). The surface gravity and the metallicity values are shown in the title of each plot.

4.3 From EW to A(Li)

With the equivalent width fitting functions FFII presented in Sect. 4.2 and defined by Eq. (7), it is in principle already possible to derive the 3D NLTE and 1D NLTE lithium abundances corresponding to a given input equivalent width by numerical inversion. Using the provided FFII procedure, one could construct a curve-of-growth for an appropriate range of abundances that maps to a range of equivalent widths covering the “observed” EW. The desired lithium abundance, or , for that particular EW may then be easily derived through interpolation in the inverted curve-of-growth EW → A(Li).

As an alternative to this numerical procedure, we have developed a third parametric fitting formula (hereafter FFIII) that represents the inverse of FFII. As such, it provides the 3D NLTE or 1D NLTE value of A(Li) for a given input EW and stellar parameters T_eff, log g, and [Fe∕H]. As an algebraic inversion of FFII is not feasible, we construct the inverse fitting function FFIII from scratch, based on the raw data (EW as a function of A(Li)).

The EW grid we start with has already been interpolated to the nominal temperature scale T_nom and completed by extrapolation where data are missing, as described in Sect. 4.2. In principle this grid would be ready for the fitting, but in this case we encounter the additional difficulty that the data in the T_nom – EW plane are not given on a rectangular grid, as shown in Fig. 7 (upper panel) for the 3D case log g = 3.5, [Fe∕H] = 0.0. While fitting data on an irregular grid is not impossible, it becomes easier if the parameter space is rectangular. Then it is possible to use the standard procedures to obtain a good fit of the data with a polynomial function of relatively low degree. However, it is not feasible to simply remap the data points to such a rectangular grid, as this would require an unacceptable amount of extrapolation. We also note that, from a physical point of view, a rectangular T_nom – EW grid is not adapted to the systematic temperature dependence of the line strength (see Fig. 7 upper panel).

In view of theses difficulties, we introduce the concept of a “scaled equivalent width” w. As a first step, we define two auxiliary fitting curves indicated in the lower panel of Fig. 7 as p(z) and q(z). These two quadratic curves approximate the EW values as a function of z = (T_nom∕1000 − 5) for the extreme abundances of our grid, A(Li) = 1.0 (p, blue curve) and A(Li) = 3.0 (q, red curve), respectively. The scaled equivalent width is then defined by the following temperature-dependent conversion: (8)

where p(z) = p₀ + p₁z + p₂z² and q(z) = q₀ + q₁z + q₂z² are the aforementioned boundary fitting curves. We note that the quality of linear fits was found to be unsatisfactory.

In principle, the coefficients p_m, q_n depend on both log g and [Fe∕H]. However, the fact that the EW values are only very weakly dependent on surface gravity allows us to consider two “global” parabolas to scale all the equivalent width values for a given metallicity. The coefficients p_m ([Fe∕H]) and q_n ([Fe∕H]) were obtained by averaging the respective individual coefficients for the three surface gravities. This approximation has the advantage of decreasing the number of numerical coefficients defining our inverse fitting function FFIII. The averaged coefficients are provided in Table B.5 as a function of [Fe∕H] for both the 3D NLTE and the 1D NLTE case.

The EW scaling largely removes the systematic trend with T_eff, as shown in the lower panel of Fig. 7, where we now plot the distribution of the data points in the z − w plane. By construction, 0 ≤ w ≤ 1, but the grid is still not perfectly rectangular. However, a rectangular grid can be easily created now by interpolation, without the need for anyextrapolation.

For this purpose, we selected 21 equidistant values between w = 0 and 1 with a step of Δw = 0.05 that will constitute, together with the temperature scale z, the rectangular parameter space on which our lithium abundances are distributed. The 21 A(Li) values that correspond to the fixed w grid are found by piecewise cubic, monotonic interpolation, separately for each z, log g, and [Fe∕H].

A this point the rectangular grid of lithium abundances can be fitted in the rectangular z − w parameter space. Also in this case, we choose to fit the surface of the abundance values by a two-dimensional polynomial function of maximum degree 4, following exactly the same method as for the abundance corrections (FFI, Sect. 4.1) and the equivalent widths (FFII, Sect. 4.2), again assuming a linear gravity dependence of the coefficients.

The final analytic expression can be written as: (9)

with w = (log₁₀(EW) − p(z))∕(q(z) − p(z)), p(z) = p₀ + p₁z + p₂z², q(z) = q₀ + q₁z + q₂z², and x = log g − 3.5, z = (T_nom∕1000 − 5).

The coefficients a_ijk are presented as a function of [Fe∕H] in Appendix B, in Tables B.6 and B.7 for the 3D NLTE and the 1D NLTE case, respectively, while the coefficients p_m and q_n describing the transformation from EW to w are provided in Table B.5.

Employing Eq. (9), optionally using the provided Python script that in addition performs interpolation in [Fe∕H], the user obtains directly the 3D NLTE or the 1D NLTE lithium abundance for given stellar parameters and the measured EW, avoiding intermediate steps such as the aforementioned inversion that is required if one uses the forward fitting function FFII for this purpose.

We have checked that our inverse fitting function FFIII works very well. The fitting function error is evaluated by taking the equivalent widths EW(T_eff, log g, [Fe∕H], A(Li)) on the grid points of our full 4D parameter space as input and comparing the output of FFIII, A(Li)_FFIII, with the true (input) lithium abundance A(Li)_inp. As shown in Fig. B.5, we find a maximum fitting function error of A(Li)_FFIII − A(Li)_inp = −0.033 dex, both for FFII_3DNLTE and FFII_1DNLTE, occurring at [Fe∕H] = − 0.5, log g = 4.0, T_eff = 5000 K, A(Li) = 3.0, where the equivalent with of the lithium doublet is ≳360 mÅ. Apart from these extreme conditions, FFIII performs much better, with a rms mean fitting error of ≈ 0.006 dex in both cases. The interpolation error is typically confined in the range |δ_FFIII|≲0.02 (see Fig. B.6).

Fig. 7 Distribution of the 3D NLTE data points for log g = 3.5, [Fe∕H] = 0.0 in the Teff – EW plane (upper panel), and in the Teff – w plane (lower panel) where the systematic trend with temperature is eliminated by normalizing the EWs with the help of the parabolic boundary curves p(z) and q(z) shown in the upper panel (see text for details).

Fig. 8 Comparison between the Li 670.8 nm curves-of-growth computed with the fitting function of Sbordone et al. (2010) for [Fe∕H] = − 2.0 (colored continuous lines) and with our function FFII (black solid lines) for the 3D NLTE (left) and 1D NLTE (right) cases, respectively, and for three values of surface gravity (log g = 3.5, 4.0, and 4.5). The vertical colored dashed lines indicate the upper limits of the lithium abundance range for which the relation of Sbordone et al. (2010) is well defined (beyond these limits the fitting function is evaluated in extrapolation). Different colors denote the different effective temperatures at which the relations are evaluated (see legend in upper left panel). The black dashed line represents the upper abundance limit up to which our fitting function is defined (A(Li) = 3.0 for all temperatures).

Fig. 9 Same as Fig. 8 but for [Fe∕H] = − 3.0, log g = 4.0 and 4.5. In this case,the upper limits to which the fitting function by Sbordone et al. (2010) is safely applicable coincides with A(Li) = 2.5 for all temperatures (vertical dashed orange line). Our fitting function is valid up to A(Li) = 3.0 (vertical dashed black line).

5 Comparison with other authors

In this section, we compare the results of our fitting functions with some of the existing work in the literature that is focused on lithium abundance corrections and parametric equations for deriving NLTE lithium abundances from EW measurements. The availability of such independent analysis performed on a common range of stellar parameters is used to check the consistency of our results and to discuss potential differences in the adopted methods.

5.1 A(Li) to EW: comparison with Sbordone et al. (2010)

The work presented in the previous sections has some important features in common with the work published by Sbordone et al. (2010). Both investigations are based on 3D CO5BOLD model atmospheres and the related 1D LHD models, and use the same codes to compute the NLTE departure coefficients (NLTE3D) and the line formation of the lithium doublet (Linfor3D). The immediate advantage of this comparison is that systematic differences in the results due to different input model atmospheres can be excluded.

In Figs. 8 and 9, we show a comparison of the results, in terms of CoGs obtained by applying both Sbordone et al. (2010)’s fitting function in the form EW(T_eff, log g,[Fe∕H],A(Li)) and our FFII(T_eff, log g, A(Li); [Fe∕H]) (Eq. (7)) to the lithium abundance range A(Li) = [1.0, 3.3]. Although our functions are valid over a wider range of metallicity, we limit the comparison to the stellar parameters that both works have in common, namely T_eff = 5500, 5900, 6300, 6500 K, log g = 3.5, 4.0, 4.5 and [Fe∕H] = − 2.0 (Fig. 8) and − 3.0 (Fig. 9).

It is worth noting that the lithium abundance range, specifically the upper boundary below which the function of Sbordone et al. is applicable, depends on the stellar parameters. This is shown in Figs. 8 and 9 by means of colored vertical dashed lines. Beyond these upper limits, their equation is not trustworthy, since it operates in complete extrapolation. On the other hand, our fitting function FFII is valid across the whole abundance range from A(Li) = [1.0 to 3.0], independently of the stellar parameters.

If we consider the abundance range where both fitting functions are well defined, we can conclude that the results are in a very good agreement. Only very subtle differences in the derived CoGs can be appreciated towards hotter models and at larger lithium abundances. This could be related to the different lithium model atom or to the choice of Sbordone et al. (2010) of using a single set of NLTE departure coefficients, computed for each model atmosphere assuming A(Li) = 2.2. This approximation can only be considered reasonable as long as the lithium line is rather weak. In the case of a stronger lithium line (larger A(Li)), the NLTE departure coefficients are no longer independent of the lithium abundance, thus introducing additional non-linearities in the partly saturated region of the curve-of-growth. In the present work, we therefore compute consistent departure coefficients for each individual A(Li).

We have verified that, on average, the departure coefficients tend to be driven closer to unity (i.e., the level populations are driven towards LTE) as the lines become optically thick. This implies that the Sbordone at al. approximation produces synthetic lithium lines that tend to be too weak for A(Li) >2.2. In case of the 3D model with T_eff = 5500 K, log g = 4.0, [Fe∕H] = − 2.0, and A(Li) = 3.0, for instance, we find EW = 226 mÅ when using the departure coefficients computed for A(Li) = 2.20, and EW = 263 mÅ when using the consistent departure coefficients computed for A(Li) = 3.00, thus Δ log(EW) ≈0.07 dex. While going in the right direction, this effect is too small to explain the discrepancies of the order of 0.3 dex seen in Fig. 8 at A(Li) = 3.0.

We conclude that these large differences are mainly due to extrapolation errors arising when the Sbordone et al. (2010) formula is applied outside its range of validity. Such errors can now be avoided by using the extended fitting function FFII which correctly tracks the impact of NLTE effects across the abundance grid for each given model atmosphere. If A(Li) falls in the range where both approximations are valid, FFII and Sbordone et al. (2010) are in very close agreement, as expected. Since the “inverse” fitting function of Sbordone et al., EW to A(Li), is the algebraic inversion of their “forward” fitting function, a comparison with our FFIII would not allow for any further conclusions.

5.2 1D NLTE lithium abundance corrections in the literature

In the context of 1D NLTE lithium abundance corrections, there are a number of other authors who contributed to this topic. We refer in this section to the work by Carlsson et al. (1994), Pavlenko & Magazzu (1996), Takeda & Kawanomoto (2005) and more recently Lind et al. (2009). They all provided a series of 1D NLTE lithium abundance corrections for the lithium resonance line at 670.8 nm, and in some cases also for the subordinate line at 610.3 nm, for a diversity of stellar parameters and lithium abundances. After briefly discussing their setups and techniques, we show a comparison of the 1D NLTE lithium abundance corrections between our results and those of the following authors.

• Carlsson et al. (1994): their atmosphere grid contains one-dimensional LTE plane-parallel MARCS models (Gustafsson et al. 1975) with T_eff between 4500 and 7500 K, log g = 0.0, 2.0, 4.0 and 4.5 and [Fe∕H] = −3.0, −2.0, −1.0, and 0.0. The NLTE departure coefficients are computed adopting a Li model atom with 21 levels and 20 bound–bound transitions, allowing them to compute 1D NLTE lithium abundance corrections for both the resonance and the subordinate lithium lines. For the latter, however, they do not provide the values of the corrections since they are smaller and show significantly less variation across the grid than for the resonance line. The 1D NLTE corrections for the Li I λ670.8 nm line are given, for all the models of their grid, through a set of coefficients for a fourth-order polynomial fit to the values of their corrections as a function of A(Li). The assumed microturbulence is not specified.

• Pavlenko & Magazzu (1996): the authors compute 1D NLTE and 1D LTE synthetic lithium line profiles and CoGs for three lithium lines: the resonance doublet (λ670.8 nm) and two subordinate lines (λ610.3 nm and λ812.6 nm). They provide 1D NLTE and 1D LTE equivalent widths in tabular form for a grid of Kurucz (1993) ATLAS9 model atmospheres for solar metallicity G-M dwarfs and subgiant stars with 3500 ≤ T_eff ≤ 6000 K and 3.0 ≤ log g ≤ 4.5 and for wide range of lithium abundances A(Li). Presumably, a fixed microturbulence parameter ξ_micro = 2.0 km s⁻¹ is adopted (details lacking). For the coolest models of their grid (3500 ≤ T_eff ≤ 4000 K), they find a very weak dependence of NLTE line profiles (and equivalent widths) on the surface gravity. In LTE, on the other hand, this dependence is quite strong, emphasizing the importance of using NLTE results in order to minimize the unrealistic effects of errors in log g on lithium abundance determinations. For comparing these results to our corrections, we converted the LTE and NLTE equivalent widths of Pavlenko & Magazzu (1996) to 1D NLTE abundance corrections (as a function of the 1D LTE abundance) by interpolation in their CoGs.

• Takeda & Kawanomoto (2005): they employ an extensive grid of 100 Kurucz (1993) ATLAS9 model atmospheres to determine 1D NLTE lithium abundances from the lithium resonance doublet line for a sample of 160 F-K dwarfs and subgiant stars of the Galactic disc. The range of stellar parameters considered is quite large, covering 5000 ≤ T_eff ≤ 7000 K, 3.0 ≤ log g ≤ 5.0, and − 1.0 ≤[Fe∕H] ≤ +0.4. A microturbulence parameter of ξ_micro = 1.5km s⁻¹ has been assumed across the entire grid of models. They provide grids of theoretical equivalent widths and the corresponding NLTE corrections for A(Li) = 0.5, 1.0, 1.5, 2.0, 2.5, 3.0 and 3.5.

• Lind et al. (2009): this work represents the most comprehensive 1D NLTE calculations of lithium abundance and corrections in late-type stars. By using a 1D MARCS (Gustafsson et al. 2008) model atmospheres grid spanning T_eff = [4000…8000] K, log g = [1.0…5.0] and [Fe∕H] = [0.0… − 3.0], they computed synthetic line profiles for the lithium lines at λ 670.8 nm and λ 610.3 nm for A(Li) = [−0.30… + 4.20]. For models with log g ≥ 3.0 the spectra were computed adopting a microturbulence parameter ξ_micro = 1.0 and 2.0km s⁻¹. For models with log g ≤ 3.0 they adopted ξ_micro = 2.0 and 5.0km s⁻¹. The LTE and NLTE equivalent widths were obtained by numerical integration over the line profiles, from which the two curves-of-growth can be constructed. The NLTE abundance corrections for each abundance were defined as the difference between the NLTE and the LTE lithium abundance that corresponds to the same equivalent width. This corresponds to the standard definition of NLTE abundance corrections and is consistent with the definition adopted in the present work (cf. Sects. 3.1 and 3.2).

The model atom used for the calculation of the NLTE departure coefficients includes the same energy levels for neutral lithium as Carlsson et al. (1994), accounting for the hyperfine structure of ⁷Li (six components for the λ670.8 nm line, three for the λ610.3 nm line). The corrections are available either by means of tables or via web interface³ that, similarly to our fitting function FFI_1DNLTE, allows a quick computation of the corrections as a function of the chosen input stellar parameters (T_eff, log g, [Fe∕H]), lithium abundance A(Li) and, in their case, microturbulence parameter ξ_micro.

In Figs. 10–14, we show the comparison between our results for the 1DNLTE A(Li) corrections (by applying FFI_1DNLTE) and those of the other authors described previously, for metallicities [Fe∕H] = 0.0, −0.5, −1.0, −2.0 and − 3.0. Each figure consists in 4 × 3 panels where different effective temperatures and gravities are shown in different rows and columns, respectively. Our results are represented by black lines, whereas the colored lines show the corrections according to the other authors, as explained by the legend of each panel. The full comparison across our grid of stellar parameters could only be done with Lind et al. (2009). The other authors cover limited or different ranges of metallicity, temperatures, and gravities and for this reason are missing in some of the comparison plots.

Our 1D NLTE corrections are based on 1D LHD spectra computed with the microturbulence parameter ξ_micro as listed in Table A.1. Since the microturbulence is an input parameter of the IDL routine provided by Lind et al. (2009)⁴, we can ensure that the comparison if fully consistent regarding the adopted value of ξ_micro, which is especially relevant for high lithium abundances where the lines are strong and potentially affected by small-scale velocity fields in the stellar photosphere. For the comparison with the other authors, full consistency regarding the choice of ξ_micro is not enforced.

Figure 10 shows the 1D NLTE corrections for solar metallicity. This plot provides a good overview of the general trend of the corrections, since such trends are qualitatively maintained also for the other [Fe∕H]. For the coolest models (T_eff=5000 and 5500 K with log g = 3.5 and 4.0) the curves drop steeply when the lithium resonance line saturates at large A(Li) values. On the other hand, for weaker lines (low A(Li)) the corrections become almost independent on the lithium abundance and assume a constant value. For T_eff =5000 K, our fitting function yields a positive 1D NLTE correction of + 0.18 dex (black line in the upper left plot of Fig. 10). The location of this knee is moving towards larger lithium abundance as the effective temperature increases, together with a flattening of the overall trend. At high metallicities ([Fe∕H] = 0.0, −0.5), the corrections are seen to be almost independent of log g.

The global trend of our 1D NLTE lithium abundance corrections is found to be in good agreement with the other works for every metallicity. Although some scatter can be appreciated especially at T_eff =5000 K, our results fit well within the other curves, indicating general consistency. A systematic shift of ~0.05 dex between our results and the others samples of corrections can be noticed from T_eff =6000 K onwards and for metallicity 0.0 (Fig. 10), − 0.5 (Fig. 11) and − 1.0 (Fig. 12). The origin of this offset could be traced back to the fact that different authors choose different model atmospheres. A different choice of the model atmospheres can indeed significantly affect the derived abundance corrections. For example, Pavlenko & Magazzu (1996) and Takeda & Kawanomoto (2005) employed the same ATLAS9 models, and their results correspond quite well with each other (Fig. 10, orange and green lines). If we compare their results with Lind et al. (2009) who adopt MARCS model atmospheres (blue line in Fig. 10), we can observe some offset in the corrections. Additionally, even different versions of the same type of models can lead to significant differences in the resulting abundance corrections. This is the case for Carlsson et al. (1994) (red line in Fig. 10) and Lind et al. (2009) who both use MARCS models in which the newer version employed by latter authors includes more line opacity. As explained by Lind et al. (2009), this causes a steeper temperature gradient in the upper part of the photosphere which, in turn, leads to a higher over-ionization. For this reason, the use of newer MARCS model atmospheres introduces a systematic shift in the overall trend of the abundance corrections (towards more positive, less negative corrections with respect to the results obtained with the old MARCS models), the size of which depends, however, on the stellar parameters.

Fig. 10 Comparison between the results of our fitting function FFI for the 1D NLTE lithium abundance corrections and those published by other authors for [Fe∕H] = 0.0, log g = 3.5, 4.0, and 4.5, Teff = 5000, 5500, 6000, and 6500 K. The black lines labeled as M19 show our results, whereas the colored lines show the corrections of other authors: C94 (Carlsson et al. 1994, red), P&M96 (Pavlenko & Magazzu 1996, orange), T&K05 (Takeda & Kawanomoto 2005, green), and L09 (Lind et al. 2009, blue).

Fig. 11 Same as Fig. 10 but for [Fe∕H] = −0.5. For this metallicity, a comparison is only possible with the results of Lind et al. 2009.

Fig. 12 Same as Fig. 10 but for [Fe∕H] = −1.0.

Fig. 13 Same as Fig. 10 but for [Fe∕H] = −2.0.

Fig. 14 Same as Fig. 10 but for [Fe∕H] = −3.0.

5.3 Impact of different model atmospheres

The effect of using MARCS models instead of 1D LHD models has been investigated by Harutyunyan et al. (2018). They adopted the same 1DLHD model atmospheres as in the present work and experienced a similar systematic offset in their NLTE corrections when comparing their results with Lind et al. (2009).

To understand the reason for this offset, they recomputed the 1D corrections by using one of the MARCS models of Lind et al. (2009) (Gustafsson et al. 2008) with parameters T_eff = 6500 K, log g = 4.0, [Fe∕H] = − 1.0. Assuming a lithium abundance of A(Li) = 2.0, they found that the resulting correction was 0.055 dex lower with respect to the value obtained by means of the corresponding 1D LHD model. Applying the same MARCS–LHD offset of − 0.055 dex over the whole range of A(Li) between 1.0 and 3.0, their modified corrections were found to be in almost perfect agreement with the work by Lind et al. (2009, see Fig. 15).

The result of this experiment suggests that the main responsible for the differences between our 1D NLTE abundance corrections and those of the other authors is the different choice of the input model atmospheres. In addition, other factors such as a different lithium model atom and the treatment of the relevant physical processes in the calculation of the NLTE departure coefficients (e.g., collisions with neutral hydrogen, charge transfer reactions) may also play a role.

6 Applications and discussion

6.1 NLTE corrections along Spite plateau

Applying our fitting functions FFI_3DNLTE and FFI_1DNLTE, we can study the influence of the different line formation models on the distribution of stars along the Spite plateau. Forthis purpose, we define a sample Spite plateau stars with stellar parameters in the range 5500 K ≤T_eff ≤6500 K (Δ = 20 K), 3.5 ≤log g ≤4.5 (Δ = 0.5), − 3.0 ≤[Fe∕H] ≤ − 1.5 (Δ = 0.025), distributed on a 3D rectangular grid. We assume that each of the 9333 stars of the sample has the same lithium abundance of A(Li) = 2.2, such that, by construction, they lie on a perfectly flat, infinitely thin Spite plateau when analyzed in 3D NLTE.

When instead analyzed in 1D LTE, the stars are distributed in a lithium abundance range 2.167 ≤ A(Li) ≤ 2.268 (see Fig. 16, star symbols). The largest negative deviation from the plateau of − 0.033 dex is found for stars with T_eff = 6500 K, log g = 3.5, [Fe∕H] = − 1.5, while the most positive deviation of + 0.068 dex is found for stars with T_eff = 5500 K, log g = 4.5, [Fe∕H] = − 2.0. The mean LTE Li abundance of the sample is unchanged, A(Li) = 2.2006, the standard deviation is as small as σ = 0.0162 dex. On the other hand, the analysis of the sample with 1D NLTE line formation is characterized by a systematic shift of the Li abundance by − 0.03 dex to A(Li) = 2.17 (see distribution of filled dots in Fig. 16). The total spread of the abundances is σ = 0.0134 dex, similar to the 1D LTE case.

We conclude from this exercise that the error introduced by analyzing the Spite plateau with 1D LTE models instead of employing the full 3D NLTE machinery is small, at most ± 0.05 dex in extreme cases. Our results also indicate that the 1D NLTE analysis may introduce a small systematic offset of ΔA(Li) ≈−0.03 dex. We see no evidence that 3D and / or NLTE effects might be responsible for the “meltdown” of the Spite plateau at metallicities below [Fe∕H] ≈ − 3.0, although we cannot demonstrate this explicitly since, strictly speaking, our fitting functions are not applicable to such low metallicities.

Fig. 15 Comparison between our 1D NLTE A(Li) corrections with the other authors for the model with Teff=6500 K, log g =4.0 and [Fe∕H] = −1.0. Our original corrections (M19, black solid line) have been shifted down (M19, black dots) by 0.055 as found by Harutyunyan et al. (2018) to account for the difference between LHD and MARCS model atmospheres. The version of the models adopted by the different authors are indicated in the legend.

Fig. 16 Lithium abundance distribution for a sample of stars with stellar parameters in the range 5500 K ≤ Teff ≤ 6500 K, 3.5 ≤ logg ≤ 4.5, − 3.0 ≤[Fe∕H] ≤ − 1.5, plotted as a function of effective temperature (top) and metallicity (bottom). By construction, all stars lie on a perfectly flat, infinitely thin Spite plateau at A(Li) = 2.2 when analyzed with 3D NLTE model spectra (dashed line). The corresponding 1D LTE (stars) and 1D NLTE (dots) abundances are derived by applying fitting functions FFI3DNLTE and FFI1DNLTE, respectively.See legend for further details.

From turn-off to the bottom of the RGB

The dependence of the corrections on T_eff and log g along an isochrone at metallicity [Fe∕H] ≈ − 2.2 is illustrated for three stars of globular cluster NGC 6397 studied by Nordlander et al. (2012), representing the turn-off point (TOP), the subgiant branch (SGB), and the base of the red giant branch (bRGB), respectively. The stellar parameters of these objects and their Li abundance derived in 1D NLTE, taken from Table 3 (temperature scale “100s”) of Nordlander et al. (2012), are reproduced in Cols. (1) – (5) of Table 1. With the help of our fitting functions FFI_1DNLTE and FFI_3DNLTE, we can easily estimate the 1D LTE and 3D NLTE lithium abundance, respectively, for the given 1D NLTE value. The results are collected in Cols. (6)–(7) of Table 1. We find that the 1D LTE abundance coincides with the 3D NLTE result for the TOP and (marginally) the SGB, but with the 1D NLTE result for the bRGB. Remarkably, the 3D NLTE − 1D NLTE abundance difference is + 0.035 dex for all three stars, such that all abundance gradients with T_eff remain perfectly unchanged. The small offset is not affecting the conclusions of Nordlander et al. (2012) whatsoever.

6.2 Microturbulence dependence

In this work, microturbulence is assumed to be a function of the fundamental stellar parameters according to Eq. (1). For this reason, ξ_micro is not considered as an independent parameter of our fitting functions.

The microturbulence dependence given by Eq. (1) is based on CO5BOLD 3D hydrodynamical simulations of the small-scale velocity field in stellar atmospheres and gives a reasonable approximation to empirical microturbulence trends measured inreals stars, as demonstrated by Dutra-Ferreira et al. (2016). It might nevertheless happen that the ξ_micro parameter measured in individual objects would exhibit major deviations from the adopted relation. In this case, our fitting functions FFII (A(Li) → EW) and FFIII (EW →A(Li)) would become unreliable for stronger lines whose equivalent width is susceptible to the influence of microturbulence. However, the abundance corrections are much less sensitive to ξ_micro (see, e.g., Lind et al. 2009), and the results obtained from fitting function FFI should still be valid, as long as the measured ξ_micro is used for the determination of the 1D LTE lithium abundance A(Li).

6.3 Blend lines

The fact that the fitting functions elaborated in this work are constructed for the pure, unblended Li feature does not imply that they cannot be applied to stars where some blend lines of other elements interfere with the lithium doublet. In this case, the EW value of the lithium resonance doublet cannot be measured directly from an observed spectrum but must be deduced by some technique todisentangle the contributions of the interfering blend lines from that of the pure lithium spectrum. This is not an uncommon procedure and there are different methods to achieve such EW extraction. For example, the solar spectrum is characterized by a weak lithium feature immersed in several blend lines. Nevertheless, Brault & Mueller (1975) succeeded in subtracting the spectrum of the blend lines (mainly the laboratory spectrum of the CN molecular band) from that of the Sun, thus isolating the lithium feature. Another technique to derive the pure Li I EW from a blended spectrum is to use the spectrum of a reference star with similar stellar parameters, but being essentially lithium-free, as a template which is subtracted from the original blended spectrum to retrieve the de-blended pure lithium line profile (e.g., Strassmeier et al. 2015; Flores Soriano et al. 2015). A necessary condition for this procedure to be valid is that all lines must be sufficiently weak, lying on the linear part of the curve-of-growth.

In general, whenever a LTE lithium abundance can be estimated from an observed spectrum, the 3D NLTE (or 1D NLTE) correction calculated with our fitting functions can be applied.

6.4 Lithium isotopic ratio

The present investigation is based on synthetic spectra that were computed for a fixed (meteoritic) lithium isotopic ratio of ⁶Li∕⁷Li = 8.2%, assuming that the derived abundance corrections are insensitive to this choice. For [Fe∕H] = 0.0 and − 2.0, we have checked the validity of this assumption by recomputing the full set of 1D LHD spectra with ⁶Li∕⁷Li = 0.0, which is probably a more realistic assumption since ⁶Li is expected to be completely destroyed during the pre-main sequence evolution of the stars in the parameter range under consideration. We find that the corrections computed with ⁶Li∕⁷Li = 0.0 start to deviate significantly from the standard corrections (⁶Li∕⁷Li = 8.2%) only for the strongest lines (EW ≳250 mÅ) due to saturation effects that lead to slightly more negative corrections. In the most extreme case, the difference amounts to − 0.038 dex (T_eff = 5000 K, log g = 3.5, [Fe∕H] = − 2.0, A(Li) = 3.0).

The fitting functions FFI, FFII, FFIII were designed to fit the data computed with ⁶Li∕⁷Li = 8.2%. However, we found that they can also be used to obtain the corrections valid for ⁶Li∕⁷Li = 0.0 and the corresponding conversion between EW and A(Li) by applying an empirical shift of the input A(Li) (FFI, FFII) or the input EW (FFIII) as follows: (10)

where quantities with subscript 0 indicate the results adjusted for zero ⁶Li contribution. With this modification, the fitting functions reproduce the 1D data computed with ⁶Li∕⁷Li = 0.0 with almost the same accuracy as the original fitting functions reproduce the data computed with ⁶Li∕⁷Li = 8.2%. We assume that the same manipulation works also for the intermediate metallicities as well as for the 3D case.

7 Conclusions

We extensively investigated the lithium line formation in solar-type FGK stars of different metallicity from a purely theoretical perspective. From the CIFIST grid of 3D CO5BOLD model atmospheres, we utilized a large sub-set of models and calculated 3D NLTE and 1D NLTE lithium abundance corrections from synthetic curves-of-growth, spanning a wide range of stellar parameters and lithium abundances. This allowed us to study the magnitude of such 3D and NLTE effects across the grid. The derived corrections can be applied to improve the accuracy of 1D LTE abundance determinations for stellar metallicities in the range − 3.0 ≤[Fe∕H] ≤ 0.0.

We found these corrections to be rather small (but still significant) for weaker lines (low and higher T_eff), where they are of the order of ± 0.05 dex. Their magnitude increases substantially for cooler temperatures and lower metallicities, especially for high lithium abundances () where the corrections can reach the − 0.30 dex level. It is also worth noting that the abundance corrections become systematically more negative towards the more metal-poor stars (with a possible reversal of the trend near [Fe/H] = − 2.0) for which reliable lithium abundances are of particular interest.

Since 3D model atmosphere simulations and NLTE line formation calculations are time consuming and require considerablecomputational resources, we first condensed the results of our study into a set of tables that provide the 3D NLTE abundance corrections as a function of the stellar parameters and the 1D LTE lithium abundance A^* (Li) to which the correction shall be applied. We then translated the tabulated corrections into analytical functions that evaluate the 3D NLTE or 1D NLTE corrections at arbitrary T_eff, log g, A^* (Li), separately for a fixed grid of metallicities [Fe∕H]. In a similar way, we also provide analytical fitting functions for directly converting a given lithium abundance into an equivalent width or vice versa, a given equivalent width into a lithium abundance (3D NLTE or 1D NLTE). The coefficients of these multi-dimensional polynomial functions are tabulated in the appendix.

For further convenience, we have developed Python scripts⁵ that perform the evaluation of all six fittingfunctions for given T_eff, log g, [Fe∕H], and ⁶Li∕⁷Li, using the polynomial coefficients tabulated asa function of [Fe∕H], applying the corrections for isotopic ratio if needed, and subsequently interpolating to the desired [Fe∕H].

By means of the fitting functions developed in this work, the results of complex 3D and NLTE calculations are readily accessible and quickly applicable to large samples of stars. In this way, our work contributes to improving the accuracy of spectroscopic lithium abundance determinations across a wide range of stellar metallicities, and may thus help to address the open questions related to the lithium content in metal-poor and solar-like stellar atmospheres.

Our results confirm that for metal-poor stars on the Spite plateau, 3D and NLTE effects are by far too small to provide a solution for the cosmological lithium problem, the persistent disagreement between the stellar plateau at A(Li) ≈ 2.2 and the primordial lithium abundance of A(Li) ≈ 2.5 inferred from WMAP and Planck measurements of the cosmic microwave background. At the same time, the metallicity dependence of the corrections offers no explanation for the Li “meltdown” observed below [Fe/H] ≈ − 3.0.

It has recently been claimed that stellar turbulent diffusion models are able to explain the lithium depletion of the Spite plateau, although their inclusion in stellar evolution models is not exempted from some uncertainties. To better constrain the underlying physics, extremely precise lithium abundance measurements along isochrones of metal-poor globular clusters are necessary. Using our fitting functions for the lithium abundance corrections, it is easy to show that differential 3D and NLTE effects between stars at the turn-off point and the base of the red giant branch are less than 0.04 dex.

Acknowledgements

We thank the State of Brandenburg (MWFK) and the German Federal Ministry of Education and Research (BMBF) for their continuous funding of basic solar-stellar research at the Leibniz Institute for Astrophysics Potsdam (AIP). A.M. is grateful for financial support from the Leibniz Graduate School for Quantitative Spectroscopy in Astrophysics, a joint project of AIP and the Institute of Physics and Astronomy of the University of Potsdam. The critical comments of an anonymous referee helped to improve the manuscript significantly.

Appendix A Model atmospheres used

In Table A.1 we compile the main properties of the full grid of models used in this work. The name of the 3D model (Col. 1) is followed by the number of opacity bins (Col. 2) and the number of representative snapshots evaluated for the respective model (Col. 3). Columns (4)–(6) list the fundamental stellar parameters (T_eff, log g and [Fe∕H]) which are identical for 3D and 1D LHD models. The effective temperature reported in Col. (4) is obtained by averaging the emergent radiative flux over the sub-sample of the selected N_snap snapshots of each CO5BOLD simulation. Finally, Col. (7) shows the microturbulence parameter ξ_micro which is computedby means of Eq. (1) and is adopted for the spectral line profile synthesis from the 1D LHD models.

Appendix B Coefficients

In this appendix we tabulate the numerical coefficients for the fitting functions provided in this work. In all of the following tables, the first column contains the name of each numerical coefficient and in the next five columns their numerical value for each one of the five metallicities of our grid. We provide Python scripts for evaluating all six fitting functions using the coefficients tabulated in this appendix online⁶.

B.1 Fitting functions FFI

Tables B.1 and B.2 list the coefficients c_ijk defined by Eq. (6) for FFI_3DNLTE and FFI_1DNLTE, respectively.Figure B.1 gives an overview of the fitting function errors at the nodes of the grid for all metallicities, while Fig. B.2 illustrates the interpolation errors evaluated on a refined temperature grid of 1D LHD models for [Fe/H] = 0.0 and − 2.0.

B.2 Fitting functions FFII

Tables B.3 and B.4 contain the coefficients w_ijk defined by Eq. (7) for the equivalent width fitting functions FFII_3DNLTE and FFII_1DNLTE, respectively.Figure B.3 gives an overview of the fitting function errors at the nodes of the grid for all metallicities, while Fig. B.4 illustrates the interpolation errors evaluated on a refined temperature grid of 1D LHD models for [Fe/H] = 0.0 and − 2.0.

B.3 Fitting functions FFIII

In the context of the “inverse” fitting function FFIII, the coefficients p_m and q_m needed for the normalization of the input equivalent width according to Eq. (8), are listed in Table B.5. The coefficients a_ijk for the related fitting function (Eq. (9)) are presented in Table B.6 for the 3D case and in Table B.7 for the 1D case. Figure B.5 gives an overview of the fitting function errors at the nodes of the grid for all metallicities, while Fig. B.6 illustrates the interpolation errors evaluated on a refined temperature grid of 1D LHD models for [Fe/H] = 0.0 and − 2.0.

Fig. B.1 Performance of the Li abundance correction fitting functions FFI3DNLTE (top) and FFI1DNLTE (bottom). For each metallicity, a total of 165 points indicates the difference (FFI− input correction) as a function of an index that runs over A(Li), Teff, and log g. Large (left group), intermediate (middle group), and small dots (right group) correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right (5000 K: red, 5500 K: black, 5900 K: gray, 6300 K: black, 6500 K: blue). Each effective temperature is represented by 10 dots indicating the results for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

Fig. B.2 Partial view of Fig. B.1 (bottom), showing for metallicities [Fe/H] = 0.0 and −2.0 the difference (FFI1DNLTE− input correction) as a function of spectrum number, now for a refined temperature grid with Δ Teff = 100 K. For each metallicity, large, intermediate, and small dots correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right. Color coding for the temperatures of the original grid is the same as in Fig. B.1; in addition, green dots represent the results for the intermediate 1D LHD test models. As before, each effective temperature is represented by 10 dots for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

Fig. B.3 Performance of the Li equivalent width fitting functions FFII3DNLTE (top) and FFII1DNLTE (bottom). For each metallicity, a total of 165 points indicates the relative difference (FFII− input EW) as a function of an index that runs over A(Li), Teff, and log g. Large (left group), intermediate (middle group), and small dots (right group) correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right (5000 K: red, 5500 K: black, 5900 K: gray, 6300 K: black, 6500 K: blue). Each effective temperature is represented by 10 dots indicating the results for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

Fig. B.4 Partial view of Fig. B.3 (bottom), showing for metallicities [Fe/H] = 0.0 and −2.0 the difference (FFII1DNLTE− input EW) as a function of spectrum number, now for a refined temperature grid with ΔTeff = 100 K. For each metallicity, large, intermediate, and small dots correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right. Color coding for the temperatures of the original grid is the same as in Fig. B.3; in addition, green dots represent the results for the intermediate 1D LHD test models. As before, each effective temperature is represented by 10 dots for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

Fig. B.5 Performance of the NLTE Li abundance fitting functions FFIII3DNLTE (top) and FFIII1DNLTE (bottom). For each metallicity, a total of 165 points indicates the difference (FFIII− input abundance) as a function of an index that runs over A(Li), Teff, and log g. Large (left group), intermediate (middle group), and small dots (right group) correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right (5000 K: red, 5500 K: black, 5900 K: gray, 6300 K: black, 6500 K: blue). Each effective temperature is represented by 10 dots indicating the results for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

Fig. B.6 Partial view of Fig. B.5 (bottom), showing for metallicities [Fe/H] = 0.0 and −2.0 the difference (FFIII1DNLTE− input abundance) as a function of spectrum number, now for a refined temperature grid with Δ Teff = 100 K. For each metallicity, large, intermediate, and small dots correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right. Color coding for the temperatures of the original grid is the same as in Fig. B.5; in addition, green dots represent the results for the intermediate 1D LHD test models. As before, each effective temperature is represented by 10 dots for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

References

All Tables

All Figures

Fig. 1 Overview of the grid of CO5BOLD 3D model atmospheres adopted for the computation of the fitting functions. The green filled squares mark the models that have been used in this work (five different metallicities, see legend at upper left), the grey squares indicate existing models that are part of the CIFIST grid but have not been considered in the present context. The blue squares highlight the missing models for which abundance corrections were derived by extrapolation (see Sect. 3.2). The solar 3D CO5BOLD model is indicated by the Sun symbol.

In the text

Fig. 2 Two different definitions of abundance corrections for the example of a strong Li line on the non-linear part of the CoG. The correction is indicated by the black horizontal arrow starting on the 3D NLTE CoG (black) at A(Li) = 2.8 and ending at a higher abundance given by the intersection with the 1D LTE CoG (blue). Conversely, the correction Δ3DNLTE-1DLTE is given by the blue horizontal arrow starting on the 1D LTE CoG (blue) at the same abundance A(Li) = 2.8 and extending to the 3D NLTE CoG (black).

In the text

Fig. 3 Overview of the Li abundance corrections Δ3DNLTE-1DLTE (top) and Δ1DNLTE-1DLTE (bottom) as evaluated on the complete 4D parameter grid. For each metallicity, a total of 165 points indicates the corrections as a function of an arbitrary spectrum number that runs over A(Li), Teff, and log g. Large (left group), intermediate (middle group), and small dots (right group) correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right (5000 K: red, 5500 K: black, 5900 K: gray, 6300 K: black, 6500 K: blue). Each effective temperature is represented by 11 dots indicating the results for 1.0 ≤ A(Li) ≤ 3.0.

In the text

	Fig. 4 Contour plots of the 3D NLTE abundance corrections in the Teff–A(Li)1DLTE plane. Each panel shows the corrections computed with FFI3DNLTE (Eq. (6)) for all combinations of surface gravity and metallicity as given in the title of each plot.
In the text

	Fig. 5 Contour plots of the 1D NLTE abundance corrections in the Teff–A(Li)1DLTE plane. Each panel shows the corrections computed with FFI1DNLTE (Eq. (6)) for all combinations of surface gravity and metallicity as given in the title of each plot.
In the text

	Fig. 6 Comparison between the input data (filled dots, color-coded according to effective temperatures as shown in the legend of the upper left panel) and the resulting fitting polynomial FFII3DNLTE (Eq. (7), red continuous line). The surface gravity and the metallicity values are shown in the title of each plot.
In the text

	Fig. 7 Distribution of the 3D NLTE data points for log g = 3.5, [Fe∕H] = 0.0 in the Teff – EW plane (upper panel), and in the Teff – w plane (lower panel) where the systematic trend with temperature is eliminated by normalizing the EWs with the help of the parabolic boundary curves p(z) and q(z) shown in the upper panel (see text for details).
In the text

Fig. 8 Comparison between the Li 670.8 nm curves-of-growth computed with the fitting function of Sbordone et al. (2010) for [Fe∕H] = − 2.0 (colored continuous lines) and with our function FFII (black solid lines) for the 3D NLTE (left) and 1D NLTE (right) cases, respectively, and for three values of surface gravity (log g = 3.5, 4.0, and 4.5). The vertical colored dashed lines indicate the upper limits of the lithium abundance range for which the relation of Sbordone et al. (2010) is well defined (beyond these limits the fitting function is evaluated in extrapolation). Different colors denote the different effective temperatures at which the relations are evaluated (see legend in upper left panel). The black dashed line represents the upper abundance limit up to which our fitting function is defined (A(Li) = 3.0 for all temperatures).

In the text

	Fig. 9 Same as Fig. 8 but for [Fe∕H] = − 3.0, log g = 4.0 and 4.5. In this case,the upper limits to which the fitting function by Sbordone et al. (2010) is safely applicable coincides with A(Li) = 2.5 for all temperatures (vertical dashed orange line). Our fitting function is valid up to A(Li) = 3.0 (vertical dashed black line).
In the text

Fig. 10 Comparison between the results of our fitting function FFI for the 1D NLTE lithium abundance corrections and those published by other authors for [Fe∕H] = 0.0, log g = 3.5, 4.0, and 4.5, Teff = 5000, 5500, 6000, and 6500 K. The black lines labeled as M19 show our results, whereas the colored lines show the corrections of other authors: C94 (Carlsson et al. 1994, red), P&M96 (Pavlenko & Magazzu 1996, orange), T&K05 (Takeda & Kawanomoto 2005, green), and L09 (Lind et al. 2009, blue).

In the text

	Fig. 11 Same as Fig. 10 but for [Fe∕H] = −0.5. For this metallicity, a comparison is only possible with the results of Lind et al. 2009.
In the text

	Fig. 12 Same as Fig. 10 but for [Fe∕H] = −1.0.
In the text

	Fig. 13 Same as Fig. 10 but for [Fe∕H] = −2.0.
In the text

	Fig. 14 Same as Fig. 10 but for [Fe∕H] = −3.0.
In the text

Fig. 15 Comparison between our 1D NLTE A(Li) corrections with the other authors for the model with Teff=6500 K, log g =4.0 and [Fe∕H] = −1.0. Our original corrections (M19, black solid line) have been shifted down (M19, black dots) by 0.055 as found by Harutyunyan et al. (2018) to account for the difference between LHD and MARCS model atmospheres. The version of the models adopted by the different authors are indicated in the legend.

In the text

Fig. 16 Lithium abundance distribution for a sample of stars with stellar parameters in the range 5500 K ≤ Teff ≤ 6500 K, 3.5 ≤ logg ≤ 4.5, − 3.0 ≤[Fe∕H] ≤ − 1.5, plotted as a function of effective temperature (top) and metallicity (bottom). By construction, all stars lie on a perfectly flat, infinitely thin Spite plateau at A(Li) = 2.2 when analyzed with 3D NLTE model spectra (dashed line). The corresponding 1D LTE (stars) and 1D NLTE (dots) abundances are derived by applying fitting functions FFI3DNLTE and FFI1DNLTE, respectively.See legend for further details.

In the text

Fig. B.1 Performance of the Li abundance correction fitting functions FFI3DNLTE (top) and FFI1DNLTE (bottom). For each metallicity, a total of 165 points indicates the difference (FFI− input correction) as a function of an index that runs over A(Li), Teff, and log g. Large (left group), intermediate (middle group), and small dots (right group) correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right (5000 K: red, 5500 K: black, 5900 K: gray, 6300 K: black, 6500 K: blue). Each effective temperature is represented by 10 dots indicating the results for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

In the text

Fig. B.2 Partial view of Fig. B.1 (bottom), showing for metallicities [Fe/H] = 0.0 and −2.0 the difference (FFI1DNLTE− input correction) as a function of spectrum number, now for a refined temperature grid with Δ Teff = 100 K. For each metallicity, large, intermediate, and small dots correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right. Color coding for the temperatures of the original grid is the same as in Fig. B.1; in addition, green dots represent the results for the intermediate 1D LHD test models. As before, each effective temperature is represented by 10 dots for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

In the text

Fig. B.3 Performance of the Li equivalent width fitting functions FFII3DNLTE (top) and FFII1DNLTE (bottom). For each metallicity, a total of 165 points indicates the relative difference (FFII− input EW) as a function of an index that runs over A(Li), Teff, and log g. Large (left group), intermediate (middle group), and small dots (right group) correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right (5000 K: red, 5500 K: black, 5900 K: gray, 6300 K: black, 6500 K: blue). Each effective temperature is represented by 10 dots indicating the results for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

In the text

Fig. B.4 Partial view of Fig. B.3 (bottom), showing for metallicities [Fe/H] = 0.0 and −2.0 the difference (FFII1DNLTE− input EW) as a function of spectrum number, now for a refined temperature grid with ΔTeff = 100 K. For each metallicity, large, intermediate, and small dots correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right. Color coding for the temperatures of the original grid is the same as in Fig. B.3; in addition, green dots represent the results for the intermediate 1D LHD test models. As before, each effective temperature is represented by 10 dots for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

In the text

Fig. B.5 Performance of the NLTE Li abundance fitting functions FFIII3DNLTE (top) and FFIII1DNLTE (bottom). For each metallicity, a total of 165 points indicates the difference (FFIII− input abundance) as a function of an index that runs over A(Li), Teff, and log g. Large (left group), intermediate (middle group), and small dots (right group) correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right (5000 K: red, 5500 K: black, 5900 K: gray, 6300 K: black, 6500 K: blue). Each effective temperature is represented by 10 dots indicating the results for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

In the text

Fig. B.6 Partial view of Fig. B.5 (bottom), showing for metallicities [Fe/H] = 0.0 and −2.0 the difference (FFIII1DNLTE− input abundance) as a function of spectrum number, now for a refined temperature grid with Δ Teff = 100 K. For each metallicity, large, intermediate, and small dots correspond to log g = 3.5, 4.0, and 4.5, respectively. For each gravity, Teff increases from left to right. Color coding for the temperatures of the original grid is the same as in Fig. B.5; in addition, green dots represent the results for the intermediate 1D LHD test models. As before, each effective temperature is represented by 10 dots for 1.0 ≤ A(Li) ≤ 2.8 plus a square for A(Li) = 3.0.

In the text