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

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.


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é et al. 2019).Lithium is the heaviest element formed in Big Bang nucleosynthesis (BBN); therefore, stellar Li abundances (A(Li) 1 ) can be used to determine the amount 12 and [X/Y] ≡ (A(X) −A(Y) ) − (A(X) −A(Y) ) ⊙ , with N X representing the number density of element "X".
of Li initially produced in BBN (Cyburt et al. 2016).Old, mainsequence 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 et al. 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 et al. 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 et al. 1996;Chabrier et al. 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 et al. 2012(Mucciarelli et al. , 2014(Mucciarelli et al. , 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 et al. 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;Bensby & Lind 2018;Aguilera-Gómez et al. 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 et al. 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 et al. 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 utilise 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.
GALactic Archaeology with HERMES (GALAH; De Silva et al. 2015) is a large spectroscopic survey conducted with the HER-MES spectrograph (Sheinis et al. 2015) using the 2dF fibre positioning system (Lewis et al. 2002) at the 3.9-metre 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 nonlocal 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 radiationhydrodynamics 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 summarise our findings.With 669 570 spectra containing 581 149 unique stars, this is the largest set of 3D NLTE Li abundances published to date.

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  = 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 (T eff , log(g), and [Fe/H]) from GALAH DR3 for this work.GALAH DR3 derived 1D NLTE T 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 T 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.

ANALYSIS
We derive the Li abundance based on direct synthetic line profile fits in a three step process: (i) 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.
(ii) 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.
(iii) 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 crossmatch, the dataset 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).

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.Each blending line is modelled by one Gaussian, given by: where A is the amplitude of the line, FWHM is the full width half max of the line, v rad is the radial velocity,  0 is the line center, and  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 database (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.

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.
First, we fit the broad region, using a single value for the blending line width FWHM and radial velocity v 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 v rad from the broad region are applied, and we now treat the continuum normalisation 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 v 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 v rad for all lines in the spectrum.The radial velocity v 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 Table 1.The wavelengths of lines used for fitting.The Li line center 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.calibration.The width of the line, FWHM, is assumed to be the same for all lines because with the typical 10 km s −1 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 v 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 v 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 v 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 v rad in stars with higher rotation due to increased blending.

Special cases
We deviate from the standard aforementioned fitting method for poorly constrained stars and stars outside of the breidablik grid (see Section 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 cutoff 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 v 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 inter-polated 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.

Equivalent width 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  2 likelihood using the nested sampling Monte Carlo algorithm, MLFriends (Buchner 2016(Buchner , 2019)), from ultranest (Buchner 2021).
ultranest is able to sample posterior distributions and return converged chains free from typical tuning parameters such as burnin 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).
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% 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%, 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.

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.minimizeadding a small error region around these parameters given by the error formula (Norris et al. 2001): where  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 normalisation parameter, we sample the posterior again over a larger region.If this re-sampling still fails to capture the full posterior, then we flag this star.Flags are discussed further in Section 3.6.

Non-detections
We determine whether Li is detected using a 2 sigma cut on the EW: EW < 2 l EW . (3) For stars where the Li line is not detected, we report the abundance corresponding to 2 l EW as the upper limit.By this definition of nondetection, 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  l 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 l EW .

Lithium abundance measurements
We estimate the Li abundance, A(Li), from the measured EW using the 3D NLTE package breidablik, taking the T 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 T eff ( T eff ) and error in EW ( EW ). A(Li) due to  T eff is provided by calculating Li abundance when perturbing T eff by  T 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 ( l EW ), and the upper  A(Li) is the Li abundance corresponding to the upper error in EW ( u 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.

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Å centered 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 normalisation 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 normalisation 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 normalise 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 +2 −5 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 normalisation factor of 0.995, which places the measured continuum higher than the true continuum by 0.5%.

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 crossmatched stars span 4280 ≤ T eff ≤ 6826 K, 1.31 ≤ log(g) ≤ 4.58, and −2.42 ≤ [Fe/H] ≤ 0.57.We crossmatch to Hypatia as this catalogue has the largest set of crossmatched 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 crossmatched 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 T 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.

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 ≤ T 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 ≤ T 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 Section 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 normalised stellar parameters falling within a sphere of radius √ 3 × 0.5 2 of the closest regularised grid point.In 3D hydrodynamic model atmospheres, T eff is an output rather than an input, so the true T eff differ from the nominal T eff , causing the grid to be irregular in T eff .To prevent holes in the grid, we use the nominal T eff values to determine whether a star is in the grid or not.To normalise the stellar parameters, we divide the stellar parameter by their respective step sizes, 500 K for T eff , 0.5 dex for log(g), and 1 for [Fe/H].
It is possible for the measured T eff to be within the grid, and the  T eff to perturb that T eff to be outside of the grid, i.e.T eff is in the grid but one or both of T eff ±  T eff is outside of the grid.Where T eff +  T eff is outside of the grid, we report  A(Li) from T eff −  T eff ; and vice versa.If both T eff ±  T eff fall outside of the grid, then  A(Li) due to  T 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.

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 utilising 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 utilising A(Li), we recommend flag_ALi < 2, removing 41939 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.
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).

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 ∼600k 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 minutes to 2 hours.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 coadded 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.We use T eff , log(g), and [Fe/H] from GALAH DR3 to re-derive 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(Mucciarelli et al. , 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).
We report the measured Li EW for all stars regardless of detection or non-detection.Fig. 6 shows the variation of Li EW across T 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 T eff ≈ 6000 K and log(g) ≈ 4.5 extending up the subgiant branch to T eff ≈ 5700 K and log(g) ≈ 4. In addition, there is a region to the left with low Li EW at T 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 et al. 1990;Deliyannis & Pinsonneault 1997), and mass loss (Schramm et al. 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 T 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 red giant branch stars.However, the median Li EW in the giant branch is low, indicating that neither population has significant Li retained.
We probe how the mean error in Li EW varies across T eff and log(g) in Fig. 7.There is a sudden increase in Li EW error at T eff above 6500 K, with some extremely high Li EW error bins scattered all throughout this high T 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.
We show how S/N affects detections and non-detections in Fig. 8.The left panel shows giant stars, using T 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 nondetections 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 cutoff 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 ≤ T 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  In general, the error increases the hotter the star, with a jump in mean error for some stars hotter than T eff ≈ 6500 K, mainly related to the rotational velocity.There is no strong correlation between log(g) and the error in Li EW.
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% of stars with no Li in our data set.The non-detections mode is centered at zero, indicating that we are separating detections and non-detections for Li-dip stars.

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 linelist 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 finetuned approach that emphasises 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 T 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 T eff .Our abundances can be higher than the GALAH DR3 abundances, this mainly occurs for giant stars at T 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%, resulting in a positive correction.The right panel shows the distribution of A(Li) for dwarfs and giants, where 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.
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 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.
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 T 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 cutoff in GALAH DR3 errors at  A(Li) ≈ 0.15 dex because stars with larger error are non-detections.This cutoff is causing the curvature at high GALAH DR3 errors.A similar cutoff is not present for our errors as we do not incorporate an error cutoff for detections.The secondary track where our errors are higher than GALAH DR3 errors are caused by large  T eff values, driving up our errors.Fig. 11 compares detections between our data set and GALAH DR3.This work detects Li in 50048 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 cutoff compared to GALAH DR3.In addition to these rapid rotators, GALAH DR3 also detects more stars with Li at T eff ≈ 4100 K, likely due to identified overdensities in the GALAH DR3 parameter space (Buder et al. 2021).

CONCLUSION
In this paper we publish 3D NLTE Li abundances and equivalent widths 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 T 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 nondetections.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 T 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 e.g.WEAVE (Dalton et al. 2016  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 cutoff in GALAH DR3 errors that is not present in our errors.The gray dashed line show where the two errors are equal.

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.
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)

Figure 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.

Figure 2 .
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 LiEW 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.

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

Figure 4 .
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 inWang et al. (2020).Our abundances are roughly in agreement with Hypatia but are generally lower even when including a 3D NLTE correction.The gray dashed line is where the abundances are equal.

Figure 5 .Figure 6 .Figure 7 .
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.

Figure 8 .
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.

Figure 9 .Figure 10 .
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 T eff and 3D NLTE A(Li), truncated at 0.3 dex.The mean A(Li) is negative with some positive A(Li) at T eff ≈ 4500 K.The increase in 3D NLTE A(Li) as T 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 centered at ∼-0.3 whilst giants mostly contribute to the distribution centered at ∼0.0.A(Li) is not only due to 3D effects, but is also caused by a different analysis method.

Figure A1 .
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: T 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.

Figure B1 .
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 .
Figure B2.Example fit of a poorly constrained star, where there are less than 3 lines above the noise in the broad region, labeled like Fig. 1.

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

Figure C1 .
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% of the pdf only intercepts the pdf once, so the upper error is reported as 34% above the MAP (green).

Figure D1 .
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.

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

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.

Table A2 .
with arbitrary values of T eff , log(g), [Fe/H], and A(Li) is also recorded.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.

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.