3D NLTE Lithium abundances for late-type stars in GALAH DR3

ABSTRACT

Lithium’s susceptibility to burning in stellar interiors makes it an invaluable tracer for delineating the evolutionary pathways of stars, offering insights into the processes governing their development. Observationally, the complex Li production and depletion mechanisms in stars manifest themselves as Li plateaus, and as Li-enhanced and Li-depleted regions of the HR diagram. The Li-dip represents a narrow range in effective temperature close to the main-sequence turn-off, where stars have slightly super-solar masses and strongly depleted Li. To study the modification of Li through stellar evolution, we measure 3D non-local thermodynamic equilibrium (NLTE) Li abundance for 581 149 stars released in GALAH DR3. We describe a novel method that fits the observed spectra using a combination of 3D NLTE Li line profiles with blending metal-line strength that are optimized on a star-by-star basis. Furthermore, realistic errors are determined by a Monte Carlo nested sampling algorithm which samples the posterior distribution of the fitted spectral parameters. The method is validated by recovering parameters from a synthetic spectrum and comparing to 26 stars in the Hypatia catalogue. We find 228 613 Li detections, and 352 536 Li upper limits. Our abundance measurements are generally lower than GALAH DR3, with a mean difference of 0.23 dex. For the first time, we trace the evolution of Li-dip stars beyond the main sequence turn-off and up the subgiant branch. This is the first 3D NLTE analysis of Li applied to a large spectroscopic survey, and opens up a new era of precision analysis of abundances for large surveys.

stars: abundances – techniques: spectroscopic – stars: late-type

1 INTRODUCTION

The chemical history of the Milky Way is encapsulated by its stars, where FGK-type stars act as fossils due to their long lifetime and convective surface which, to first approximation, retains the abundances they were born with (Jofré, Heiter & Soubiran 2019). Lithium is the heaviest element formed in Big Bang nucleosynthesis (BBN); therefore, stellar Li abundances (A(Li) ¹) can be used to determine the amount of Li initially produced in BBN (Cyburt et al. 2016). Old, main-sequence turn-off (MSTO) dwarf stars have been found to exhibit the same Li abundance over a range of metallicities, which is a phenomenon known as the Spite plateau (Spite & Spite 1982). The Spite plateau abundance of Li, at A(Li) ≈2.2 dex (Bonifacio & Molaro 1997; Ryan, Norris & Beers 1999; Asplund et al. 2006; Meléndez et al. 2010; Sbordone et al. 2010) is a factor of 3 lower than the 2.75 dex predicted by BBN (Pitrou et al. 2018); this inconsistency is referred to as the cosmological Li problem (Fields 2011). This difference is likely due to gravitational settling and turbulent mixing in stars depleting and depositing Li below the convective surface (Richard, Michaud & Richer 2005; Korn et al. 2007). Through the cosmological Li problem, Li links BBN and stellar evolution.

Lithium is used to probe stellar evolution because it is fragile, burning at temperatures above 2.5 MK (Basri, Marcy & Graham 1996; Chabrier, Baraffe & Plez 1996; Bildsten et al. 1997; Ushomirsky et al. 1998). Li depletes when a star undergoes the first dredge-up due to the deepening of the convective envelope evolving from the main sequence onto the red giant branch (RGB). Li is further sharply depleted at the RGB bump (Lind et al. 2009b). Due to these two depletion events, metal-poor giants also form a plateau (Mucciarelli, Salaris & Bonifacio 2012; Mucciarelli et al. 2014, 2022), similar to the Spite plateau. The Li-dip is a depletion of Li for MSTO stars within a narrow temperature range of 6400–6850 K (Boesgaard & Tripicco 1986; Gao et al. 2020), first observed in the Hyades cluster (Wallerstein, Herbig & Conti 1965), later observed in other clusters (Balachandran 1995; Burkhart & Coupry 2000) and field stars (Randich et al. 1999; Chen et al. 2001; Lambert & Reddy 2004; Ramírez et al. 2012; Aguilera-Gómez, Ramírez & Chanamé 2018; Bensby & Lind 2018), with the largest set of field stars observed from GALAH DR3 in Gao et al. (2020). Models are able to reproduce the observed depletion of Li in the Hyades cluster through atomic diffusion, rotation induced mixing, and internal gravity waves (Talon & Charbonnel 1998; Montalbán & Schatzman 2000). To understand how these Li depletion mechanisms vary with stellar properties such as age and metallicity requires accurate Li abundance measurements for a large set of stars.

Elemental abundances cannot be directly measured, instead they must be inferred from models through radiative transfer: the propagation of radiation through the stellar atmosphere (Asplund 2005; Lind & Amarsi 2024). Radiative transfer can be solved under the assumption of local thermodynamic equilibrium (LTE), which is carried out using the Saha and Boltzmann equations to calculate electron level populations and opacities. Non-LTE (NLTE) line formation is the relaxation of LTE and influences the line opacity by including the radiation field in calculating level populations. This can change the atmospheric depth at which an absorption line is formed, producing a more realistic line shape compared to LTE line formation. Measured abundances therefore differ in NLTE compared to LTE depending on the element, stellar parameters, and line strength. For Li, the difference between measured LTE and NLTE abundance can be up to 0.4 dex (Lind, Asplund & Barklem 2009a). Model stellar atmospheres are used in radiative transfer to simulate the medium that radiation propagates through. 1D hydrostatic model atmospheres are often used and approximate convection through mixing-length theory, micro-, and macroturbulence. These approximations are relaxed in the more realistic 3D hydrodynamic model atmospheres, which model convection from first principles (Stein & Nordlund 1998; Freytag et al. 2012). Again, measured abundances differ in 3D compared to 1D depending on the element, stellar parameters, and line strength (Collet, Asplund & Trampedach 2006; Dobrovolskas et al. 2013). For Li, the difference between measured 1D and 3D abundances can be up to 0.2 dex (estimated from data published in Wang et al. 2020).

Although 3D and NLTE effects are often discussed separately, in 3D NLTE radiative transfer these effects are coupled; meaning that the difference between measured abundance in 1D LTE compared to 3D NLTE is not a straightforward sum of the difference in LTE and NLTE with the difference in 1D and 3D (Klevas et al. 2016; Lind & Amarsi 2024). For Li, 3D and NLTE effects partially cancel in MSTO stars, resulting in corrections close to zero when compared to 1D LTE abundances. However, these effects do not cancel in giant stars, resulting in up to 0.5 dex difference in measured Li abundance (Wang et al. 2020). 3D NLTE models provide a Li line that matches observed Li lines more closely in shape compared to 1D LTE Li profiles, but were unfeasible to utilize outside of small samples of stars due to the computational cost associated with such models (Lind et al. 2013; Mott et al. 2017; Wang et al. 2022). However, 3D NLTE models are now becoming more accessible with large grids being produced and published (Ludwig et al. 2009; Magic et al. 2013). Due to the simplicity of the Li atom (see Wang et al. 2020 and the references therein), there are two such grids currently available for Li (Harutyunyan et al. 2018; Wang et al. 2020). Through these grids, it is now possible to apply 3D NLTE measurements for Li to large spectroscopic surveys.

Existing large-scale catalogues of Li in the literature include AMBRE, GALAH DR2, LAMOST, and Gaia-ESO. AMBRE measures 1D LTE Li abundances with NLTE corrections from Lind et al. (2009a) for ∼7000 stars from FEROS, HARPS, and UVES. These spectra were degraded to the lowest resolution of the three spectrographs, R =40 000 (Guiglion et al. 2016). GALAH DR2 measures 1D NLTE Li abundances using spectra from HERMES at R =28 000 for ∼340 000 stars using the data driven Cannon approach (Buder et al. 2018). LAMOST measures 1D LTE Li abundances for ∼160 000 stars at R ≈7500. Gaia-ESO measures 1D LTE Li abundances derived from curves of growth using spectra from UVES (R =47 000) and GIRAFFE (R ≈17 000) for ∼38 000 stars (Franciosini et al. 2022). All existing surveys explore a similar parameter space of FGKM-type stars, with most surveys stopping near [Fe/H] =−2.5, whereas AMBRE goes lower to [Fe/H] =−5. Notably, none of the existing Li catalogues are in 3D NLTE.

GALactic Archaeology with HERMES (GALAH; De Silva et al. 2015) is a large spectroscopic survey conducted with the HERMES spectrograph (Sheinis et al. 2015) using the 2dF fibre positioning system (Lewis et al. 2002) at the 3.9-m Anglo-Australian Telescope. Data release 3 (DR3) provided stellar parameters and abundances of 30 elements for 669 570 spectra (Buder et al. 2021). These parameters were derived through Spectroscopy Made Easy (sme; Valenti & Piskunov 1996; Piskunov & Valenti 2017) with 1D marcs model atmospheres (Gustafsson et al. 2008), using non-local thermodynamic equilibrium (NLTE) radiative transfer for 12 elements (Amarsi et al. 2020), including Li. However, this analysis does not consider 3D effects such as the impact of the radiation–hydrodynamics interaction on the atmospheric structure; nor does it consider the interaction between 3D and NLTE effects which are non-linear.

In this paper, we present 3D NLTE Li abundances based on GALAH DR3 spectra and stellar parameters, derived using 3D hydrodynamic model atmospheres under NLTE radiative transfer, with improved treatment of blends and error estimates. In Section 2, we discuss the GALAH observations and parameters used in this work. In Section 3, we present the analysis method used to derive 3D NLTE Li abundances and errors. In Section 4, we present our measured 3D NLTE Li abundances. In Section 5, we compare our new abundances to the GALAH DR3 abundances. Lastly, in Section 6 we summarize our findings. With 669 570 spectra containing 581 149 unique stars, this is the largest set of 3D NLTE Li abundances published to date.

2 OBSERVATIONS AND STELLAR PARAMETERS

We measure Li abundances from spectra observed for GALAH DR3, using 3D NLTE spectrum synthesis. GALAH has 4 CCDs each with discrete wavelength channels, covering 4713–4903 Å, 5648–5873 Å, 6478–6737 Å, and 7585–7887 Å, respectively (De Silva et al. 2015). We only consider the spectra on CCD 3 as this is where the Li 6707.814 Å line of interest falls. The resolution in CCD 3 varies with wavelength and from fibre to fibre over the range |$R = 22\, 000$|–32 000, and is typically ∼25 500 in the Li region (Kos et al. 2017).

The spectra were reduced, processed, and analysed as described by Kos et al. (2017) and Buder et al. (2021). We retain the stellar parameters (⁠|$\rm {T}_{\rm {eff}}$|⁠, log(g), and [Fe/H]) from GALAH DR3 for this work. GALAH DR3 derived 1D NLTE |$\rm {T}_{\rm {eff}}$| and [Fe/H] based on SME synthetic spectra and measured log(g) from Gaia DR2 (Lindegren et al. 2018) and Hipparcos (van Leeuwen 2007) parallaxes (see Buder et al. 2021 for details). They compared their stellar parameters to accurate photometric estimates and found very good agreement for cool stars, while temperatures were underestimated by upwards of 125 K for the warmest stars. Although the |$\rm {T}_{\rm {eff}}$| values in GALAH DR3 are not based on the most accurate 3D NLTE Balmer line analyses (see Amarsi et al. 2018), they appear to be close to the fundamental scale aside from issues with the warmest stars.

We use S/N adopted from snr_c3_iraf, representing the S/N per pixel measured in CCD3 which contains the Li line. In addition to these parameters, we also adopt the GALAH DR3 flags associated with these parameters. In total we analysed 669 570 spectra, reporting on 581 149 stars, using stacked spectra where repeated observations of the same star were taken.

3 ANALYSIS

We derive the Li abundance based on direct synthetic line profile fits in a three step process:

• Each spectrum is fitted in a two-step process, over a broad region then a narrow region. The narrow Li region is fit with a 3D NLTE breidablik (Wang et al. 2020) line profile for Li, and Gaussian line profiles that represent several blending lines.

• From this fit, we derive the equivalent width (EW) of the best-fitting Li line, along with lower and upper errors using a Monte Carlo nested sampling algorithm.

• We then evaluate the Li abundance, A(Li), at this EW using breidablik and stellar parameters from GALAH DR3 (Buder et al. 2021).

breidablik is a set of interpolation routines made available for a grid of Li line profiles calculated in 3D NLTE that were published in Wang et al. (2020). We update breidablik for this work as discussed in Appendix  A. The following subsections discuss these steps in more detail, special cases where we deviate from this process, and other technical details of the analysis.

We have chosen a traditional spectrum fitting approach for this work as opposed to label transfer: transferring Li abundances from an existing catalogue to GALAH DR3 spectra. Whilst label transfer through the use of neural networks has been fast and accurate for spectroscopic surveys on Li (Nepal et al. 2023), there is not a large enough training set of 3D NLTE Li abundances with blending lines to create a sufficiently accurate model for this task. However, given a large enough cross-match, the data set produced in this paper can become a new training set for a label transfer model of 3D NLTE Li abundances. This opens up a new method of measuring 3D NLTE Li abundances for future surveys like WEAVE (Dalton et al. 2016) and 4MOST (de Jong et al. 2019).

3.1 Spectrum fit

3.1.1 Model spectra

The Li line is heavily blended for solar metallicity stars, which are the majority of stars in GALAH. For each spectrum, we fit a 3D NLTE Li line profile and Gaussian absorption lines to blending lines in the region. The overall model spectrum is then produced by multiplying these theoretical spectra together.

We use a 3D NLTE profile for Li because it is non-Gaussian in shape when saturated. These profiles are from breidablik, which contains a grid of 3D NLTE L think i line profiles interpolated to arbitrary stellar parameters through radial basis functions (RBFs). Each blending line is modelled by one Gaussian, given by

$$\begin{eqnarray} \rm {flux} = 1-\rm {A} \exp {\left(\frac{(\lambda -\lambda _0 (1+\rm {v}_{\rm {rad}}/c))^2}{FWHM^2/4\ln (2)}\right)}, \end{eqnarray}$$

(1)

where |$\rm {A}$| is the amplitude of the line, FWHM is the full width at half-maximum of the line, |$\rm {v}_{\rm {rad}}$| is the radial velocity, λ₀ is the line centre, and c is the speed of light. The breidablik Li line profiles are broadened using a Gaussian convolution with their own separate width (FWHM_Li). This is because breidablik line profiles are the result of a detailed radiative transfer solution through a 3D hydrodynamic model atmosphere, and thus have a significant intrinsic width due to stellar atmospheric effects and quantum mechanics (including thermal, convective, and microscopic collisional broadening), whilst Gaussian absorption lines have an intrinsic width of zero.

The line lists for the blending lines were compiled from Meléndez et al. (2012) and the VALD3 data base (Ryabchikova et al. 2015), supplemented with newer molecular data; the wavelengths and original sources are provided in Table 1. There is a V line at 6708.094 Å and Ce line at 6708.099 Å; these lines have very similar strengths and are too close to be resolved, so are modelled as one line. The Fe i 6713.742 Å line is two Fe i lines merged together. There is a prominent spectral feature mainly in giant stars at ∼6709 Å that we could not identify the species for. We fit 2 Gaussians to this feature, labelled as ‘?1’ and ‘?2’ in Table 1.

Table 1.

The wavelengths of lines used for fitting. The Li line centre is taken as the weighted sum of Li components. V and Ce are at very similar strengths and wavelengths and are therefore modelled as one line.

Broad region		Narrow region
Centre	Reference	Centre	Reference
Al i 6696.085	Eriksson & Isberg (1963)	CN 6706.730	Mandell, Ge & Murray (2004)
Al i 6698.673	Eriksson & Isberg (1963)	Fe i 6707.433	Nave et al. (1994)
Fe i 6703.565	Nave et al. (1994)	CN 6707.545	Mandell et al. (2004)
Fe i 6705.101	Nave et al. (1994)	Li i 6707.814	Smith, Lambert & Nissen (1998)
Fe i 6710.317	Nave et al. (1994)	V/Ce 6708.096	Meggers, Corliss & Scribner (1975)/Palmeri et al. (2000)
Fe i 6711.819	Nave et al. (1994)	?1 6708.810	This work
Fe i 6713.095	Nave et al. (1994)	?2 6709.011	This work
Fe i 6713.742	Nave et al. (1994)
Ca i 6717.681	Risberg (1968)

3.1.2 Spectrum fit

Each spectrum is fitted in two steps using a broad (6695–6719 Å) and narrow (6706–6710 Å) region. The fits to both regions are shown in Fig. 1.

Figure 1.

Fit to an example spectrum of a solar metallicity giant. Fitting of both the broad 6695–6719 Å region (left) and the narrow 6706–6710 Å region (right) corresponding to the green shaded region in the left panel is shown. Each theoretical spectrum is coloured per element, with the vertical lines showing the central wavelength of the line. The red shaded region in the right panel is the range that the synthetic spectrum would change by if there was a 1σ change in EW. The combined theoretical spectrum is shown in blue. Spectral properties in the bottom left are from GALAH DR3.

First, we fit the broad region, using a single value for the blending line width FWHM and radial velocity |$\rm {v}_{\rm {rad}}$| based on nine Al, Fe, and Ca lines in the broad region, omitting the narrow region, shown in the left panel of Fig. 1. Then, we fit the narrow region, where the fitted FWHM and |$\rm {v}_{\rm {rad}}$| from the broad region are applied, and we now treat the continuum normalization as a variable offset. We model the blending lines with six Gaussians, whilst Li is modelled with breidablik line profiles broadened by FWHM_Li, shown in the right panel of Fig. 1. We adopt this two-step approach to constrain the |$\rm {v}_{\rm {rad}}$| and FWHM free from influence from the Li line, thus reducing the number of free parameters fitted over the narrow region. The 6706.730 CN line is not modelled perfectly by a Gaussian; however, because it overlaps so little with the Li line, this will not affect our derived Li EW significantly. Additional example fits of a warm solar metallicity dwarf with low-S/N, a poorly constrained star, and a saturated star are shown in Appendix  B.

We fit a single value of FWHM and |$\rm {v}_{\rm {rad}}$| for all lines in the spectrum. The radial velocity |$\rm {v}_{\rm {rad}}$| represents the residual velocity of the Li line after the spectrum has been shifted to the rest frame, primarily due to a velocity drift caused by errors in the wavelength calibration. The width of the line, FWHM, is assumed to be the same for all lines because with the typical 10 km s⁻¹ resolution of the HERMES instrument, the line shape is dominated by the instrumental profile together with rotational broadening rather than broadening intrinsic to the line or due to convection. The fitted FWHM and |$\rm {v}_{\rm {rad}}$| from this work are slightly different from the values reported in GALAH DR3 as we are only fitting these values over the 6695–6719 Å region as opposed to GALAH which fits over all CCDs; and we define these values differently. In particular, FWHM from this work represents the combined rotational and instrumental broadening under the assumption of a Gaussian line profile.

To prevent unconstrained fits in stars with no measurable Li but a depressed continuum, we apply constraints on both FWHM and |$\rm {v}_{\rm {rad}}$|⁠. We take the rotational broadening value measured from GALAH DR3, and assume that the instrumental resolution lies in the range R=22000–32000 (Kos et al. 2017); and sum the rotational and instrumental broadening in quadrature as limits for FWHM. A similar upper limit is placed for FWHM_Li, but due to the intrinsic broadening of breidablik line profiles, the lower limit is 0. We limit |$\rm {v}_{\rm {rad}}$| to prevent it from fitting noise in the continuum when there are no lines present. This limit is empirically set to half the upper FWHM limit because it is harder to measure |$\rm {v}_{\rm {rad}}$| in stars with higher rotation due to increased blending.

3.1.3 Special cases

We deviate from the standard aforementioned fitting method for poorly constrained stars and stars outside of the breidablik grid (see subsection 3.5 for details on the breidablik grid).

A spectrum is classified as poorly constrained if fewer than three lines can be measured in the broad region. This most often occurs when stars are metal-poor, hot, or rapidly rotating. The three line cut-off is determined empirically. If a spectrum is poorly constrained, then the Li line is only fit with a breidablik line profile assuming no contribution from other elements. The FWHM_Li and |$\rm {v}_{\rm {rad}}$| are fitted simultaneously with Li EW as there is no strong constraint on these values from the broad region.

We also deviate from the standard fitting method for stars with stellar parameters outside of the breidablik grid. These stellar parameters are too far from the breidablik grid and so the interpolated 3D NLTE Li line profile is untrustworthy. Instead, we use a Gaussian to also model the Li line, with the same constraints as the blending lines. These Li EWs are reliable except for stars with saturated Li lines.

3.2 EW measurements

The EW of the Li line is calculated via numerically integrating the best-fitting breidablik line profile. To account for measurement error introduced by uncertainties in the strengths of the blending lines, we sample the posterior distribution of the fit given by χ² likelihood using the nested sampling Monte Carlo algorithm, MLFriends (Buchner 2016, 2019), from ultranest (Buchner 2021).

ultranest is able to sample posterior distributions and return converged chains free from typical tuning parameters such as burn-in and the number of samples. This is possible because ultranest samples a posterior distribution by iteratively shrinking the parameter space towards higher likelihoods and iteratively updating the evidence. In order to capture the full posterior distribution and to distinguish between detections and non-detections, we allow Li to be in emission. This is implemented as a reflection of the absorption line with the same EW. The blending lines can only be in absorption or else the solution is poorly constrained.

To interpret the sampled posterior, we make a histogram using the sampled posterior, with bin widths given by the Stone algorithm (Stone 1984) because it produces reasonable bin widths for both normal and skewed distributions. We then smooth the histogram with a running mean of window size 3 to reduce noise. From this density histogram, we measure the Li EW and errors as shown in Fig. 2. The measured Li EW is adopted as the maximum a posteriori (MAP) estimate (the mode of the smoothed histogram) as this typically fits the observed spectra better than the mean; and the lower and upper errors are given by the highest posterior density Bayesian credible interval (Hespanhol et al. 2019). This method of analysing sampled posteriors is very similar to Hinton (2016).

Figure 2.

Example sampled posterior showing how we derive statistics. The sobject_id and S/N of this star is shown in the top right corner. The histogram is the sampled posterior binned (black bars), with the smooth density (solid black line) obtained via a running mean. The measured Li EW is adopted as the MAP (red); and the lower and upper errors are adopted as the highest posterior density Bayesian credible interval (blue), where the horizontal line used to identify the error in Li EW is also shown.

Fig. 2 shows our lower and upper errors on a well-behaved distribution. The errors are given by the highest posterior density Bayesian credible interval: the EWs with the same probability density value that encompass a 68 per cent region of the sampled posterior. To find these EWs, first, a horizontal line is drawn, then the EWs where it intersects the sampled posterior are identified, lastly the area encompassed by these EWs are measured. This horizontal line is adjusted until the area encompassed by the EWs is 68 per cent, and the corresponding EWs are then the 1σ lower and upper errors. The highest posterior density Bayesian credible interval reduces to a standard ±1σ confidence interval when the distribution is symmetric and unimodal. There are 170 stars where the sampled posterior is not fully captured due to restrictions from the prior, we discuss how the errors for these stars are derived in Appendix  C.

3.2.1 Sampled posterior probability distribution region

The posterior is sampled over a bounded region. We set this bounded region using the parameters of the best fit from scipy.optimize.minimize adding a small error region around these parameters given by the error formula (Norris, Ryan & Beers 2001)

$$\begin{eqnarray} \sigma _{\rm {EW}} = \lambda n^{1/2}_{\rm {pix}}/R (\textrm {S/N}), \end{eqnarray}$$

(2)

where R is the resolving power, S/N is the signal-to-noise ratio per pixel, and n_pix is the number of pixels in the line. The posterior is then well sampled if most of the probability distribution function is captured. We define this by considering where the mode of the sampled posterior falls. First, we take the histogram of the sampled posterior for one parameter using 100 bins, then we identify the mode of this histogram. If the mode is within 5 bins of either the lower or upper edge of the histogram, then this sampled posterior is not well captured for this parameter. This process is then repeated for all parameters. If the posterior is not well captured for any EWs or the continuum normalization parameter, we sample the posterior again over a larger region. If this resampling still fails to capture the full posterior, then we flag this star. Flags are discussed further in subsection 3.6.

3.2.2 Non-detections

We determine whether Li is detected using a 2 sigma cut on the EW

$$\begin{eqnarray} \rm {EW} \lt 2 \sigma ^{\rm l}_{\rm {EW}}. \end{eqnarray}$$

(3)

For stars where the Li line is not detected, we report the abundance corresponding to |$2 \sigma ^{\rm l}_{\rm {EW}}$| as the upper limit. By this definition of non-detection, negative EWs are also necessarily non-detections that we report an upper limit for. We note that errors in the model are not taken into account, so our errors are underestimated.

In cases where |$\sigma ^{\rm l}_{\rm {EW}}$| does not exist, reported as NaN in the data (see right panel of Fig. 2), we the Norris et al. (2001) error formula to determine whether the line is detected. In the case of a non-detection we report upper limit 2σ_EW, similar to stars with a non-NaN error, 2|$\sigma ^{\rm l}_{\rm {EW}}$|⁠.

3.3 Lithium abundance measurements

We estimate the Li abundance, A(Li), from the measured EW using the 3D NLTE package breidablik, taking the |$\rm {T}_{\rm {eff}}$|⁠, log(g), and [Fe/H] values from GALAH DR3. breidablik interpolates from a grid of 3D NLTE Li EWs and A(Li) to arbitrary stellar parameters using a feedforward neural network. There are two main sources of error for A(Li) (σ_A(Li)): error in |$\rm {T}_{\rm {eff}}$| (⁠|$\sigma _{\rm {T}_{\rm {eff}}}$|⁠) and error in EW (σ_EW). σ_A(Li) due to |$\sigma _{\rm {T}_{\rm {eff}}}$| is provided by calculating Li abundance when perturbing |$\rm {T}_{\rm {eff}}$| by |$\sigma _{\rm {T}_{\rm {eff}}}$| from GALAH DR3 and reported as the mean of the lower and upper error. σ_A(Li) due to σ_EW is provided separately as an asymmetric error, where the lower σ_A(Li) is the Li abundance corresponding to the lower error in EW (⁠|$\sigma ^{\rm l}_{\rm {EW}}$|⁠), and the upper σ_A(Li) is the Li abundance corresponding to the upper error in EW (⁠|$\sigma ^{\rm u}_{\rm {EW}}$|⁠). We do not report the error in A(Li) due to error in log(g) or error in [Fe/H], as these are much smaller than the aforementioned sources of error, typically of the order of 0.01 dex.

3.4 Validation

3.4.1 Synthetic spectra

We validate our analysis pipeline on a randomly generated synthetic spectrum, intended to be similar to a heavily blended solar spectrum. To create this synthetic spectrum, we synthesize a broad region spectrum and narrow region spectrum separately then multiply them together. The broad region spectrum is created through our template spectrum using EWs of ∼70 mÅ. In addition, we include 40 Gaussians of EW ∼20 mÅ centred at random wavelengths in the broad region excluding the narrow region to mimic the potential impact of blending lines that we do not model. The narrow region is then created through the template spectrum with EWs of ∼20 mÅ. We add a continuum normalization factor of 0.99 and noise corresponding to S/N=100. The fitted synthetic spectrum and true parameters are shown in Appendix  D.

We show each component of the synthetic spectrum along with the fits in Fig. 3. For this comparison, we remove the continuum normalization factor from the synthetic spectrum for clarity. The left panel shows that the measured EWs are higher for all lines compared to the true EWs, and this is worse where unmodelled blends are stronger, showing unmodelled blends inflate the measured EWs. Unmodelled blends also cause the measured FWHM to be higher than the true FWHM. In addition, we do not normalize the broad region spectrum, which contributes to the overestimated parameters for this example. Whilst this does affect our Li EW measurements, it does not introduce a clear bias towards higher or lower Li EWs. Unlike the broad region, the fitted narrow region EWs can either be smaller or larger than the true EWs. instead, the fitted spectrum better follows the noisy spectrum as opposed to the synthetic spectrum, as shown in the right panel, implying that the differences in measured EW for the narrow region is due to noise rather than blending lines. We fit a Li EW of |$17_{-5}^{+2}$| mÅ, which gives A(Li) =1.74 dex; compared to the true Li EW of 20 mÅ corresponding to A(Li) =1.80 dex. We note however that the true value almost falls within 1σ of our measured Li EW, despite a large portion of the offset being caused by systematic rather than random errors. We recover a continuum normalization factor of 0.995, which places the measured continuum higher than the true continuum by 0.5 per cent.

Figure 3.

The components to the synthetic spectrum and the fitted spectrum. The left panel shows the broad region without noise or continuum normalization factor, with the unmodelled blends shown. The right panel shows the narrow region without continuum normalization factor, for both the synthetic spectrum with and without noise.

3.4.2 Hypatia catalogue

We validate that our treatment of blends successfully recovers the Li line by comparing to a sample of stars from the Hypatia catalogue (Hinkel et al. 2014), shown in Fig. 4. The cross-matched stars span 4280 ≤ |$\rm {T}_{\rm {eff}}$| ≤6826 K, 1.31 ≤ log(g) ≤4.58, and −2.42 ≤ [Fe/H] ≤0.57. We cross-match to Hypatia as this catalogue has the largest set of cross-matched stars with lithium in GALAH, but note that as the Hypatia data are collected from different sources there may be significant systematics present. For these 26 cross-matched stars, we find similar but slightly lower abundances for 1D LTE A(Li) from Hypatia. we derive 3D NLTE abundance corrections from breidablik (Wang et al. 2020), and find that these remove a small amount of the offset, but is not enough to account for the full difference. We find a mean difference of 0.18 dex, with scatter 0.21 dex; applying 3D NLTE corrections, the mean difference becomes 0.14 dex whilst scatter remains unchanged. These abundances are not within 1σ agreement of each other, but it’s unclear if this is because errors are underestimated on our side or in the Hypatia catalogue. We do not find any trends in abundance difference with |$\rm {T}_{\rm {eff}}$|⁠, log(g), [Fe/H], or A(Li), indicating that our results perform equally well for all stellar parameters and are not biased by the different amounts of blending present for different stars.

Figure 4.

3D NLTE A(Li) from this work compared to A(Li) for 26 stars in common with the Hypatia (Hinkel et al. 2014) catalogue. Blue points are the 1D LTE A(Li) from Hypatia, orange points are these abundances in 3D NLTE, derived using the 3D NLTE corrections published in Wang et al. (2020). Our abundances are roughly in agreement with Hypatia but are generally lower even when including a 3D NLTE correction. The grey dashed line is where the abundances are equal.

3.5 Synthetic spectra grid

The synthetic spectra grids used in breidablik and GALAH are different, and as a result, stellar parameters derived in GALAH DR3 may be outside the breidablik grid. The GALAH DR3 analysis is based upon the 1D marcs grid (Gustafsson et al. 2008) and spans 2500 ≤ |$\rm {T}_{\rm {eff}}$| ≤8000 K, −0.5 ≤ log(g) ≤5.5 dex, and −5.0 ≤ [Fe/H] ≤1.0; whilst the breidablik grid is based on the 3D stagger grid (Magic et al. 2013) and spans 4000 ≤ |$\rm {T}_{\rm {eff}}$| ≤7000 K, 1.5 ≤ log(g) ≤5.0 dex, and −4.0 ≤ [Fe/H] ≤0.5. The Li abundance for stars outside the stellar grid are reported but flagged as untrustworthy, see subsection 3.6 for details.

For stellar parameters near the edge of the grid, it is not obvious whether this is inside or outside the grid. We consider a star to be within the breidablik grid if a star is within half the distance between two adjacent grid points to the nearest grid model. In practice, this is defined as the normalized stellar parameters falling within a sphere of radius |$\sqrt{3 \times 0.5^2}$| of the closest regularized grid point. In 3D hydrodynamic model atmospheres, |$\rm {T}_{\rm {eff}}$| is an output rather than an input, so the true |$\rm {T}_{\rm {eff}}$| differ from the nominal |$\rm {T}_{\rm {eff}}$|⁠, causing the grid to be irregular in |$\rm {T}_{\rm {eff}}$|⁠. To prevent holes in the grid, we use the nominal |$\rm {T}_{\rm {eff}}$| values to determine whether a star is in the grid or not. To normalize the stellar parameters, we divide the stellar parameter by their respective step sizes, 500 K for |$\rm {T}_{\rm {eff}}$|⁠, 0.5 dex for log(g), and 1 for [Fe/H].

It is possible for the measured |$\rm {T}_{\rm {eff}}$| to be within the grid, and the |$\sigma _{\rm {T}_{\rm {eff}}}$| to perturb that |$\rm {T}_{\rm {eff}}$| to be outside of the grid, that is, |$\rm {T}_{\rm {eff}}$| is in the grid but one or both of |$\rm {T}_{\rm {eff}}$| ±|$\sigma _{\rm {T}_{\rm {eff}}}$| is outside of the grid. Where |$\rm {T}_{\rm {eff}}$| +|$\sigma _{\rm {T}_{\rm {eff}}}$| is outside of the grid, we report σ_A(Li) from |$\rm {T}_{\rm {eff}}$| −|$\sigma _{\rm {T}_{\rm {eff}}}$|⁠; and vice versa. If both |$\rm {T}_{\rm {eff}}$| ±|$\sigma _{\rm {T}_{\rm {eff}}}$| fall outside of the grid, then σ_A(Li) due to |$\sigma _{\rm {T}_{\rm {eff}}}$| could not be evaluated and is reported as NaN. This occurs for 34 stars in our data set.

The grid covers a range in A(Li) from −0.5 dex to 4 dex, however, we still report A(Li) and σ_A(Li) outside of these limits, as extrapolation in this parameter is relatively well behaved. However, we do caution that these extrapolated abundances should be verified before use.

3.6 Flags

We use a bitmask flag (flag_ALi) to indicate details with the analysis shown in Table 2, similar to GALAH DR3. For example, if a particular star is a non-detection and the Li EW posterior is not fully captured, then flag_ALi will be 5. We recommend using flag_ALi < 4 when utilizing Li EWs, removing 219 stars, because Li EWs are well-measured if the continuum is well placed and the Li EW posterior is fully captured. When utilizing A(Li), we recommend flag_ALi < 2, removing 41 939 stars outside of the breidablik grid. This stricter flag selection is because A(Li) is sensitive to extrapolation in the breidablik grid whilst Li EW is not significantly affected.

Table 2.

Flags used in this analysis to calculate the final bitmask flag flag_ALi.

flag_ALi
1	Non-detection, upper limit reported
2	Stellar parameters fall outside the breidablik grid
4	Li EW posterior is not fully captured
8	Continuum posterior is not fully captured
16	Bad stellar parameter flag (≥ 128)

We reproduce the GALAH flags on stellar parameters, metallicity, and A(Li) that are adopted in this work in Table 3. In general, we recommend taking flag_sp_DR3 and flag_fe_h_DR3 as 0, and flag_ALi_DR3 < 2 when using the GALAH DR3 data replicated in this work. In particular, flag_ALi_DR3 = 0 for detections and flag_ALi_DR3 = 1 for upper limits. For further details on the GALAH DR3 flags, see Buder et al. (2021).

Table 3.

Flags from GALAH DR3 that are adopted from Buder et al. (2021) in this work.

flag_sp_DR3
1	Unreliable astrometric solution, Gaia RUWE > 1.4
2	Unreliable broadening
4	Low-S/N (below 10 for CCD 2)
8	Reduction issues
16	t-SNE projected emission features
32	t-SNE projected binaries
64	Binary sequence/pre-main sequence flag
128	S/N-dependent high sme χ2 (bad fit)
256	Problems with Fe
512	sme did not finish
1024	Stellar parameters fall outside of the marcs grid
flag_fe_h_DR3
1	Upper limit
2	Bad fit/large χ2
4	Uncertain measurement/saturation
16	Bad stellar parameter flag (≥ 128)
32	No measurement available
flag_ALi_DR3
1	Upper limit
2	Bad fit/large χ2
4	Uncertain measurement/saturation
8	Bad wavelength solution/rv for Li6708
16	Bad stellar parameter flag (≥ 128)
32	No measurement available

4 RESULTS

In this work, we present a comprehensive catalogue of Li abundance measurements as a supplementary catalogue for GALAH DR3. The catalogue is tabulated as described in Table 4. In total, we used ∼600 000 CPU hours to compute this catalogue, with majority of the computational time going to posterior sampling. Due to the different number of parameters fitted for poorly constrained stars compared to the standard analysis procedure, the time per star varies between 30 min and 2 h. We present both the ‘allspec’ and ‘allstar’ catalogues from GALAH DR3. ‘allspec’ contains all spectra, including repeated observations of the same star; whilst ‘allstar’ includes analysis from only the co-added spectra for stars that have been observed multiple times (refer to Buder et al. 2021 for the treatment of duplicate spectra). All analysis of results and comparison to GALAH DR3 for this work has been done with the ‘allstar’ catalogue.

Table 4.

Schema for the outputs of both the allstar and allspec catalogues for this work. The results from this work are made available through a Value Added Catalogue to GALAH DR3. We provide both results from this work and include relevant results from GALAH DR3.

Name	Data type	Units	Description
sobject_id	int	–	Spectrum ID from GALAH DR3
star_id	str	–	Star ID from GALAH DR3
fwhm_broad	float	km s−1	Width of the lines measured from broad region, NaN if poorly constrained
fwhm_Li	float	km s−1	Width of the breidablik Li line, NaN if outside grid
vbroad_DR3	float	km s−1	Rotational broadening from GALAH DR3 (no instrumental profile)
delta_rv_6708	float	km s−1	Relative |$\rm {v}_{\rm {rad}}$| compared to GALAH DR3 measured from the broad region, measured from Li region if poorly constrained
rv_DR3	float	km s−1	|$\rm {v}_{\rm {rad}}$| from GALAH DR3
EW	float	mÅ	Measured Li EW
e_EW_low	float	mÅ	Lower error in Li EW (⁠|$\sigma ^{\rm l}_{\rm {EW}}$|⁠)
e_EW_upp	float	mÅ	Upper error in Li EW (⁠|$\sigma ^{\rm u}_{\rm {EW}}$|⁠)
e_EW_norris	float	mÅ	Error using formula from Norris et al. (2001)
ALi	float	dex	3D NLTE A(Li) detections
ALi_upp_lim	float	dex	3D NLTE A(Li) upper limit corresponding to |$2\sigma ^{\rm l}_{\rm {EW}}$|
e_ALi_low	float	dex	Lower error in A(Li) associated with |$\sigma ^{\rm l}_{\rm {EW}}$|
e_ALi_upp	float	dex	Upper error in A(Li) associated with |$\sigma ^{\rm u}_{\rm {EW}}$|
e_ALi_teff	float	dex	Error in A(Li) associated with |$\sigma _{\rm {T}_{\rm {eff}}}$|
flag_ALi	int	–	Bitmask flag for this analysis (see subsection 3.6)
ALi_DR3	float	dex	1D NLTE A(Li) detections and upper limits from GALAH DR3
e_ALi_DR3	float	dex	σA(Li) from GALAH DR3
flag_ALi_DR3	int	–	Bitmask flag for A(Li) from GALAH DR3
teff_DR3	float	K	|$\rm {T}_{\rm {eff}}$| from GALAH DR3
e_teff_DR3	float	K	|$\sigma _{\rm {T}_{\rm {eff}}}$| from GALAH DR3
logg_DR3	float	log(cm s−2)	log(g) from GALAH DR3
e_logg_DR3	float	log(cm s−2)	σlog(g) from GALAH DR3
flag_sp_DR3	int	–	Bitmask flag for |$\rm {T}_{\rm {eff}}$| and log(g) from GALAH DR3
fe_h_DR3	float	dex	[Fe/H] from GALAH DR3
e_fe_h_DR3	float	dex	σ[Fe/H] from GALAH DR3
flag_fe_h_DR3	int	–	Bitmask flag for metallicity from GALAH DR3
snr_DR3	float	–	S/N of the spectra from GALAH DR3

We use |$\rm {T}_{\rm {eff}}$|⁠, log(g), and [Fe/H] from GALAH DR3 to rederive the lithium abundance. Stars that do not have stellar parameters reported in GALAH DR3 are not included in our abundance analysis. We carry over the GALAH DR3 flags and add our own flag, but we report Li EW and A(Li) where possible regardless of these flags. In total, we measure A(Li) in 669 570 spectra which contains 581 149 unique stars, of which 228 613 are detections and 352 536 are upper limits.

Our derived 3D NLTE A(Li) over [Fe/H] are shown in Fig. 5 with detections on the left panel and non-detections with A(Li) reported as upper limits on the right panel. For stars with Li detections, we see an overdensity of warm and cool dwarfs for solar metallicity stars, labelled with ovals; and we see two plateaus: the metal-poor cool dwarf plateau, also known as the Spite plateau (Spite & Spite 1982; Sbordone et al. 2010), and the metal-poor giant plateau (Mucciarelli et al. 2012; Mucciarelli et al. 2014, 2022); labelled with lines. These plateaus and overdensities are shown for illustrative purposes and not analysed further in this work. In the Spite plateau, detections and upper limits are well separated, with upper limits mostly below detections; this is in contrast to the metal-poor giant plateau, which has detections and upper limits at similar abundances. This indicates that for metal-poor dwarfs, we are distinguishing between stars with detected and non-detected Li successfully; whilst for metal-poor giants we are not. To distinguish between detections and non-detections for metal-poor giants, we need to lower the upper limits which can be achieved through higher S/N spectra. Given that detections and upper limits are very close-in abundance for giants, we find that the current GALAH data cannot be used to constrain population-synthesis modelling of Li abundance (similar to Chanamé et al. 2022).

Figure 5.

The distribution of 3D NLTE A(Li) detections (left) and upper limits (right) against [Fe/H]. There is an overdensity of warm and cool dwarfs, labelled as ovals. The cool plateau (also known as Spite plateau) and the metal-poor giant plateau are also labelled. There is also an overdensity of stars with non-detections in a similar location to the cool dwarfs with detections. A large number of metal-poor giants have upper limits at similar abundances to the metal-poor giant plateau.

We report the measured Li EW for all stars regardless of detection or non-detection. Fig. 6 shows the variation of Li EW across |$\rm {T}_{\rm {eff}}$| and log(g), with the left panel showing the mean and the right panel showing the median. Both panels show a region in the MSTO with high Li EW at |$\rm {T}_{\rm {eff}}$| ≈ 6000 K and log(g) ≈ 4.5 extending up the subgiant branch to |$\rm {T}_{\rm {eff}}$| ≈ 5700 K and log(g) ≈ 4. In addition, there is a region to the left with low Li EW at |$\rm {T}_{\rm {eff}}$| ≈ 6500 K and log(g) ≈ 4.2; the well-known Li-dip. This region extends beyond the MSTO into the sub giant branch, and is likely the evolved Li-dip population. The transport of material out of the convective layer into the core of the star depletes Li and is stronger for cooler stars with a deeper convective zone. The Li-dip forms due to the moderation of this transport which inhibits Li depletion for stars cooler than the dip. The exact mechanism which inhibits this transport is unknown (Aguilera-Gómez et al. 2018), with internal gravity waves (Charbonnel & Talon 1999; Montalbán & Schatzman 2000), diffusion (Michaud 1986; Richer & Michaud 1993), rotationally induced mixing (Pinsonneault, Kawaler & Demarque 1990; Deliyannis & Pinsonneault 1997), and mass-loss (Schramm, Steigman & Dearborn 1990) being possible explanations. The Li-dip has been seen before in GALAH data (Gao et al. 2020), but this is the first time that a population of stars that have evolved onward from the Li-dip has been seen. Li is detected again for warm dwarfs hotter than the Li-dip, however the Li EW in stars hotter than the Li-dip is lower than for stars cooler than the dip. For stars hotter than the Li-dip, the mean is sometimes negative and lower than the median indicating that Li is not detected in a small number of stars, likely rapid rotators. The Li EW is weaker in giants compared to the MSTO, however, there are some scattered bins with a higher EW in the mean which are not present in the median. This feature is due to Li rich giants, which bring up the mean significantly. The red clump is visible in the mean Li EWs at |$\rm {T}_{\rm {eff}}$| ≈ 4700 K and log(g) ≈ 2.5 dex, indicating that there are a higher proportion of Li rich red clump stars compared to Li rich RGB stars. However, the median Li EW in the giant branch is low, indicating that neither population has significant Li retained.

Figure 6.

Mean EW (left) and median EW (right) binned in |$\rm {T}_{\rm {eff}}$| and log(g) for all stars in our sample, including non-detections. Both the mean and the median show stars with a higher Li EW in the MSTO at |$\rm {T}_{\rm {eff}}$| ≈ 6000 K and log(g) ≈ 4.5, extending up the subgiant branch to |$\rm {T}_{\rm {eff}}$| ≈ 5700 K and log(g) ≈ 4. Li-dip stars are hotter with lower Li EW at |$\rm {T}_{\rm {eff}}$| ≈ 6500 K and log(g) ≈ 4.2, also extending up the subgiant branch. Li then increases in even hotter stars at |$\rm {T}_{\rm {eff}}$| ≈ 7000 K and log(g) ≈ 4, however the Li EW does not reach as high as before the Li-dip. The mean can be lower than the median for rapid rotators at |$\rm {T}_{\rm {eff}}$| ≈ 7500 K; and the mean tends to be higher than the median for giant stars, likely due to the influence of Li-rich stars.

We probe how the mean error in Li EW varies across |$\rm {T}_{\rm {eff}}$| and log(g) in Fig. 7. There is a sudden increase in Li EW error at |$\rm {T}_{\rm {eff}}$| above 6500 K, with some extremely high Li EW error bins scattered all throughout this high |$\rm {T}_{\rm {eff}}$| MSTO region. This feature is due to rapid rotators, which have increased line-blending due to high-rotational broadening, thus making the line harder to measure as the depression from the line is lower. There is no strong trend in log(g), as the evolutionary state of the star influences the amount of line-blending, we conclude that blending lines do not affect our Li EW error significantly. This is likely because we do not include model errors in our analysis.

Figure 7.

Mean error in Li EW, given by |$(\sigma ^{\rm l}_{\rm {EW}} + \sigma ^{\rm u}_{\rm {EW}})/2$|⁠. In general, the error increases the hotter the star, with a jump in mean error for some stars hotter than |$\rm {T}_{\rm {eff}}$| ≈ 6500 K, mainly related to the rotational velocity. There is no strong correlation between log(g) and the error in Li EW.

We show how S/N affects detections and non-detections in Fig. 8. The left panel shows giant stars, using |$\rm {T}_{\rm {eff}}$| < 5500 K and log(g) > 3 dex. An increase in S/N will increase the number of detections at lower EW and shift the overlap between detections and non-detections down in EW. Both these effects indicate that we are sensitive to weaker lines at higher S/N. The overlap in EW between detections and non-detections has a smaller range at higher S/N, indicating that we are better able to separate out detections and non-detections. However, the non-detections mode is slightly positive regardless of S/N, caused by some stars with a non-zero Li EW, indicating that we are not fully separating detections and non-detections for giant stars. There is a long detection tail at high Li EW comprised of Li rich giants. The smoothness of this tail indicates that there is no distinct cut-off between Li rich and Li normal giants, instead, giants exhibit the full range of Li abundances. The right panel shows stars within the Li-dip, taken as 5900 ≤ |$\rm {T}_{\rm {eff}}$| ≤ 7000 K and 3.5 ≤ log(g) ≤ 5 dex. Li-dip detections show a bimodality that is more evident at higher S/N because at higher S/N stars with less Li are detected. These high-S/N detections are likely due to false positives which appear to outnumber the true positives significantly at these low EWs. False positives are when stars with no Li are measured with Li due to the noise in the spectra. We use a 2σ detection limit, so this will occur for 2.5 per cent of stars with no Li in our data set. The non-detections mode is centred at zero, indicating that we are separating detections and non-detections for Li-dip stars.

Figure 8.

The probability density function of Li EW for giant stars (left) and Li-dip stars (right), showing both detections (solid) and non-detections (dashed). An additional S/N cut of 100 is made and shown in red. The grey dash dotted lines are at 0. Both giant and MSTO stars appear more well separated in detections and non-detections at higher S/N. However, giant stars detections are unimodal with a positive tail, whilst MSTO stars have bimodal detections.

5 DISCUSSION

We compare our 3D NLTE Li abundances to 1D NLTE Li abundances in GALAH DR3. The flags used for this comparison are flag_sp_DR3 = = 0 and flag_fe_h_DR3 = = 0. Li EW that fall outside of the breidablik grid are still accurate but these 3D NLTE Li abundances are extrapolated, so we use flag_ALi < 4 for Li EWs and flag_ALi < 2 for 3D NLTE Li abundances. For 1D NLTE abundances from GALAH DR3, we use flag_ALi_DR3 < 2.

The mean difference (δA(Li)) in A(Li) between this work and GALAH DR3 is shown in Fig. 9. The differences here are not only due to 3D effects, but also due to a difference in analysis technique. GALAH DR3 fitted synthetic spectra with sme, which uses a line list to set the blending line strengths with a scaled solar composition, and only reported A(Li) where there were more than 5 unblended Li pixels. In contrast, we have fitted spectra using a fine-tuned approach that emphasizes the impact of blending lines and include errors in blends through sampling the posterior distribution. We report A(Li) regardless of number of blended Li pixels. The left panel shows that in general our 3D NLTE abundances are lower than the GALAH DR3 abundances. This negative difference is primarily due to the difference in analysis technique and model spectra as the 1D NLTE to 3D NLTE correction for A(Li) ≈2 dex is positive for effectively the entire HR diagram (derived using data from Wang et al. 2020). The Li abundance increases as |$\rm {T}_{\rm {eff}}$| increases due to the detection limit, where a higher A(Li) is required to form a line of the same Li EW at higher |$\rm {T}_{\rm {eff}}$|⁠. Our abundances can be higher than the GALAH DR3 abundances, this mainly occurs for giant stars at |$\rm {T}_{\rm {eff}}$| ≈ 4500 K, for a range of Li abundances. Furthermore, our positive differences do show a slight trend in Li abundance, where the difference is more positive for higher Li abundance. These stars with a positive difference tend to be giant stars with a line depression of ∼20 per cent, resulting in a positive correction. The right panel shows the distribution of δA(Li) for dwarfs and giants, where we take log(g) < 3.5 dex as giants and log(g) ≥ 3.5 dex as dwarfs. On average, our abundances differ from GALAH DR3 by −0.23 dex. The all stars distribution shows a bimodality in difference: at ∼-0.3, due to both giants and dwarfs; and at ∼0.0, primarily due to giants. Similar to the left panel, this shows that on average our differences are negative, however, we do have a small number of positive differences primarily for giant stars.

Figure 9.

The difference between 3D NLTE A(Li) from this work and 1D NLTE A(Li) from GALAH DR3, defined as δA(Li). The left panel shows how δA(Li) varies as a function of |$\rm {T}_{\rm {eff}}$| and 3D NLTE A(Li), truncated at 0.3 dex. The mean δA(Li) is negative with some positive δA(Li) at |$\rm {T}_{\rm {eff}}$| ≈ 4500 K. The increase in 3D NLTE A(Li) as |$\rm {T}_{\rm {eff}}$| increases is due to the detection limit. Right panel shows δA(Li) as a histogram for giants compared to dwarfs. Here giants are stars with log(g) < 3.5 dex and dwarfs are log(g) ≥ 3.5 dex. δA(Li) forms a bimodal distribution, with giants and dwarfs both contributing to the distribution centred at ∼-0.3 whilst giants mostly contribute to the distribution centred at ∼0.0. δA(Li) is not only due to 3D effects, but is also caused by a different analysis method.

We compare our errors in A(Li) against GALAH DR3 in Fig. 10 for stars with detected A(Li). To derive the error in A(Li), 3D NLTE σ_A(Li), we take the mean error due to EW, summed in quadrature with the error in |$\rm {T}_{\rm {eff}}$|⁠. The median error in this work is similar to GALAH DR3, at σ_A(Li) ≈ 0.1 dex. However, there is a curvature in the errors where we derive higher errors for stars with both low- and high-GALAH DR3 errors. This curvature at low-GALAH DR3 errors is due to our errors having a Li EW dependence whilst GALAH DR3 errors do not. There is a sharp cut-off in GALAH DR3 errors at σ_A(Li) ≈ 0.15 dex because stars with larger error are non-detections. This cut-off is causing the curvature at high-GALAH DR3 errors. A similar cut-off is not present for our errors as we do not incorporate an error cut-off for detections. The secondary track where our errors are higher than GALAH DR3 errors are caused by large |$\sigma _{\rm {T}_{\rm {eff}}}$| values, driving up our errors.

Figure 10.

Comparison of error in A(Li) derived in this work to GALAH DR3. The errors are comparable, with the median at 0.1 dex for both data sets. There is a slight curvature in the comparison, where for low and high errors in GALAH DR3, we derive higher errors. There is a secondary track where we derive higher errors than GALAH DR3. In addition, there is a sharp cut-off in GALAH DR3 errors that is not present in our errors. The grey dashed line show where the two errors are equal.

Fig. 11 compares detections between our data set and GALAH DR3. This work detects Li in 50 048 more stars compared to DR3 over all stellar parameters, shown in the left panel. The increase in detections is in part due to differences in the detection criteria used and the difference in model spectra. The density of our detections for dwarf stars are similar to Fig. 6, indicating that we detect more stars with Li where the mean Li EW is higher. This is not true for giants, where we detect more stars with Li in the red clump than GALAH DR3, however these stars do not have a higher mean Li EW. Detections in GALAH DR3 but not this work traces a similar shape, as shown in the right panel. However, GALAH DR3 detects more stars with Li in hot MSTO stars, likely the rapid rotators where our analysis has a stricter detection cut-off compared to GALAH DR3. In addition to these rapid rotators, GALAH DR3 also detects more stars with Li at |$\rm {T}_{\rm {eff}}$| ≈ 4100 K, likely due to identified overdensities in the GALAH DR3 parameter space (Buder et al. 2021).

Figure 11.

Stars that have detections in this work but not in DR3 (left) compared to stars that have detections in DR3 but not this work (right). This work increases detections for all stellar parameters compared to DR3, however, we increase detections most in RC stars, and MSTO stars that are cooler or hotter than the Li dip. Compared to DR3, we lose detections mostly in the hot MSTO stars, where the rapid rotators are.

6 CONCLUSION

In this paper, we publish 3D NLTE Li abundances and EWs for 581 149 stars observed in GALAH DR3. For these stars, we measure Li abundances for 228 613 stars and provide upper limits for 352 536 stars. We fit the Li region of the observed spectra using model spectra produced by multiplying together a 3D NLTE Li profile and Gaussian profiles for other lines. We provide errors that consider blending lines by sampling the posterior for all our stars using these model spectra. We validate this method on synthetic spectra and by comparison to 26 stars in the Hypatia catalogue. We identify the known Li phenomena, including the Spite plateau formed by cool metal-poor dwarfs, the metal-poor giant plateau, an overdensity of warm and cool dwarfs at solar metallicity, and the Li-dip. Using our Li EWs, we find stars that evolved from the Li-dip on the main sequence turn-off into the subgiant branch. We show that the error in Li EW increases sharply for rapid rotators and the error in Li abundance has a |$\rm {T}_{\rm {eff}}$| dependence. Non-detections are defined through our Li EW errors, which are correlated with S/N. We show that higher S/N better separates detections and non-detections. Comparing our abundances with 1D NLTE GALAH DR3 Li abundances, on average we differ by 0.23 dex in measured abundance, with majority of our abundances lower than GALAH DR3. We find the detection limit, which is strongly influenced by EW, increases as |$\rm {T}_{\rm {eff}}$| increases. In addition, we find that our errors are similar to GALAH DR3, with median error of σ_A(Li) ≈ 0.1 dex for both. We detect Li in more stars than GALAH DR3, most of these stars fall in regions with higher mean Li EW; whilst GALAH DR3 detects Li in more rapid rotators compared to this work.

In future work, we plan on implementing this method for the upcoming GALAH data release 4 and using it to investigate various Li phenomena identified during our analysis. We note that the methodology presented in this paper is applicable with minimal modification to future surveys, for example WEAVE (Dalton et al. 2016) and 4MOST (de Jong et al. 2019).

ACKNOWLEDGEMENTS

We thank the anonymous referee for their many suggestions improving the quality and content of this manuscript. This work was supported by computational resources provided by the Australian Government through the National Computational Infrastructure (NCI) under the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme (project y89) and HPC-AI Talent Programme Scholarship (project hl99). This work was supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. SLM acknowledges support from the Australian Research Council through Discovery Project grant DP220102254, and from the UNSW Scientia Fellowship Programme. KL acknowledges funds from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement no. 852977) and funds from the Knut and Alice Wallenberg Foundation. We acknowledge the traditional owners of the land on which GALAH data were taken, the Gamilaraay people, and pay our respects to elders past, present, and emerging.

7 DATA AVAILABILITY

The data produced as part of this work, described in Section 4, can be downloaded as a Value Added Catalogue from the GALAH DR3 website, https://www.galah-survey.org/dr3/the_catalogues/#3d-nlte-li-abundances. The sampled posteriors can be provided upon request. The code described in Section 3 can be found at https://github.com/ellawang44/galah-li.

Footnotes

1

We use the customary abundance notation where |$\rm {A(Li)} \equiv \log (N_\mathrm{X}/N_\textrm { }H) + 12$| and [X/Y] ≡ (A(X) − A(Y)) − (A(X) − A(Y))_⊙, with N_X representing the number density of element ‘X’.

APPENDIX A: breidablik

We update the line profile interpolation routine in breidablik to use RBFs instead of Kriging, taking the implementation of RBF from Bertran de Lis et al. (2022). The methodology used for comparison is the same as presented in Wang et al. (2020).

To make the interpolation easier, we transform the flux using

$$\begin{eqnarray} f_{\rm {t}} = \log _{10}(1-f + s), \end{eqnarray}$$

(A1)

where s is a small positive constant that is treated as a hyperparmeter in the interpolation methods. For RBF, we find s = 10⁻⁵ using 5-fold cross-validation.

We compare the performance of RBF against Kriging, taking the Kriging results from Wang et al. (2020). The leave-one-out cross-validation errors, interpolation, extrapolation, training, and execution time of the code is shown in Table A1. Both the mean absolute deviation (MAD), which is less affected by outliers, and the root mean squared (rms), which is more affected by outliers, are reported. The 3D leave-one-out cross-validation errors are likely overestimates as the step sizes in this test were 500 K |$\rm {T}_{\rm {eff}}$|⁠, 0.5 dex log(g), and 0.5 dex [Fe/H], whilst in practice, interpolating on the full 3D grid will use step sizes closer to 250 K |$\rm {T}_{\rm {eff}}$|⁠, 0.25 dex log(g), and 0.25 dex [Fe/H]. To account for this, we also measure leave-one-out cross-validation errors for a 1D grid, as doubling the size of our 3D grid to do these tests is computationally infeasible. These 1D tests include 1000 interpolation points, where the data points are fully bounded by existing grid points, and 1000 extrapolation points, where the data points are not fully bounded by the existing grid. These results show that RBF is comparable to Kriging for the 3D grid and 1D interpolation, slightly better for 1D extrapolation, and is faster to train and execute.

Table A1.

Comparison of Kriging and RBF interpolation errors when measuring A(Li) from the 670.8 nm line profile. The MAD and rms error statistics for A(Li) are reported using leave-one-out cross-validation on the 3D grid, and interpolation and extrapolation tests on a 1D grid. The time required to train the model and the time required to produce the three Li lines (610.4 nm, 670.8 nm, 812.6 nm) with arbitrary values of |$\rm {T}_{\rm {eff}}$|⁠, log(g), [Fe/H], and A(Li) is also recorded.

	3D leave-one-out		1D interpolation		1D extrapolation
Method	MAD	rms	MAD	rms	MAD	rms	Training time (s)	Execution time (s)
Kriging	0.014	0.022	0.020	0.038	0.019	0.143	235	5.62
RBF	0.012	0.018	0.022	0.039	0.016	0.062	56.3	0.0273

To verify the final interpolation model and to demonstrate RBF interpolated line profiles and the usual expected errors, we measured the abundances of 4 benchmark stars. Again, we compare these results to Kriging from Wang et al. (2020). The measured abundances are are shown in Table A2, and the corresponding interpolated line profiles for the simulated A(Li) are shown in Fig. A1. For these benchmark stars, RBF outperforms Kriging slightly.

Figure A1.

Comparison between RBF and Kriging on verification models (Balder) that were not trained on.These models were tailored to specific stars, indicated by name and stellar parameters: |$\rm {T}_{\rm {eff}}$|⁠, log(g), [Fe/H], and A(Li). Note that the residuals between prediction and verification model have been magnified by a different amount in each plot.

Table A2.

A comparison between Kriging and RBF interpolation on benchmark stars. Direct synthesis in the first row is the reference A(Li), subsequent rows is the error relative to the reference. The verification models were tailored to specific stars, indicated by name, whose stellar parameters are given in Fig. A1.

Method	Sun	HD 140283	Procyon	HD 84937
Model
Direct synthesis	1.100	2.000	1.000	2.200
Kriging	0.012	0.002	−0.006	0.009
RBF	0.010	0.000	0.007	−0.004

APPENDIX B: ADDITIONAL FITS

As examples, we provide addition fits to: a warm solar metallicity dwarf with low-S/N (Fig. B1), a poorly constrained fit (Fig. B2), and a saturated fit (Fig. B3).

Figure B1.

Example fit of a warm solar metallicity dwarf with low-S/N, where the lithium line dominates over the blending lines. Labels follow Fig. 1.

Figure B2.

Example fit of a poorly constrained star, where there are less than 3 lines above the noise in the broad region. Labels follow Fig. 1.

Figure B3.

Example of a strongly saturated star, indicating an excellent fit with the synthetic spectrum. Labels follow Fig. 1.

APPENDIX C: POSTERIOR ERROR

We do not fully capture the sampled posterior for 170 stars due to restrictions from the prior, shown in Fig. C1. The highest posterior density Bayesian credible interval does not exist for these sampled posteriors. In these cases, we report 34 per cent below or above the MAP as the lower or upper error, depending on which one is defined. The corresponding undefined lower or upper error is then set to NaN.

Figure C1.

Example sampled posterior showing how we derive errors for a posterior with one side truncated. The sampled posterior is still binned into a histogram (black bars) and smoothed (solid black line), with the MAP (red) adopted as the measured Li EW. However, the horizontal line (blue) encompassing 68 per cent of the Probability Distribution Function (PDF) only intercepts the PDF once, so the upper error is reported as 34 per cent above the MAP (green).

APPENDIX D: SYNTHETIC SPECTRUM VALIDATION

We validate our method on synthetic spectrum. Fig. D1 shows the synthetic spectrum that we test on and our corresponding best fit.

Figure D1.

Example fit to a synthetic spectrum with added noise. This spectrum is intended to be similar to a heavily blended solar spectrum. Labels follow Fig. 1.

The input parameters and the fitted results are shown in Table D1.

Table D1.

The true and measured parameters for validation on a synthetic spectrum. The measured parameters are close to the true parameters, with the broad region measured parameters usually larger than the true parameters.

Parameter	Units	True	Measured
Broad region
Al i 6696.085	mÅ	75	83
Al i 6698.673	mÅ	108	122
Fe i 6703.565	mÅ	82	112
Fe i 6705.101	mÅ	77	79
Fe i 6710.317	mÅ	87	93
Fe i 6711.819	mÅ	74	96
Fe i 6713.095	mÅ	104	110
Fe i 6713.742	mÅ	110	134
Ca i 6717.681	mÅ	119	130
FWHM	km s−1	17.13	17.22
|$\rm {v}_{\rm {rad}}$|	km s−1	3.52	3.38
Narrow region
CN 6706.730	mÅ	7.2	9.1
Fe i 6707.433	mÅ	13	5.3
CN 6707.545	mÅ	19	29
Li i 6707.814	mÅ	20	17
V/Ce 6708.096	mÅ	8.4	13
?1 6708.810	mÅ	1.3	4.1
?2 6709.011	mÅ	1.6	0.3
FWHMLi	km s−1	13.68	10.09