© ESO, 2016

1. Introduction

The chemical composition of extremely metal-poor stars offers unique insight into the early evolution of the Galaxy and the properties of the first stars, the so-called Population intial mass function III. No Population III star has ever been observed, and their initial mass function (IMF) is debated due to the sensitivity of the models to the physics of accretion and fragmentation.

Several groups have run hydrodynamical, cosmological simulations coupled with primordial gas chemistry, where dark matter halos of 10⁶M_⊙ collapse gravitationally to form the first stars at redshifts z ≈ 20−30 (e.g., Gao et al. 2007; Bromm 2013). Abel et al. (2002) simulated such a minihalo, in which a self-gravitating dense core of ~100 M_⊙ cooled by molecular hydrogen, formed a single protostellar core. Depending on the accretion scenario, they argued that within a few thousand years after the simulations ended, the star would reach a mass of at least 30 M_⊙, and possibly more than 100 M_⊙ within 1 Myr. No other stars would form due to the immense radiation field of the young hot star, which would rapidly photodissociate all molecular hydrogen in the entire minihalo (pausing further star formation), and later disrupt it as the star explodes in a supernova after a few million years.

Refined simulations by Greif et al. (2011) highlight the importance of fragmentation in the accretion disk, leading to the formation of clusters of ~10 protostars per minihalo. Their simulations indicated low protostellar masses (0.1−10 M_⊙) in a flat (top-heavy) distribution. The final mass of most of these objects would likely be higher, but the accretion phase of nearly half of the objects ended already at this stage when dynamical interactions ejected them from the star forming region, meaning long-lived low-mass Population III stars could very well still exist today.

Hirano et al. (2014) ran 2D axisymmetric radiation hydrodynamical stellar structure calculations, simulating the evolution of the central protostar in a self-collapsing gas cloud until the point where UV feedback shut off accretion and the star reached the main sequence. Their simulations resulted in a top-heavy stellar mass distribution with roughly equal numbers in the ranges 10−100 and 100−1000 M_⊙. Lower-mass stars were not produced, as their simulations did not consider accretion disk fragmentation. Susa et al. (2014) reached similar conclusions from 3D simulations, with a mass spectrum covering 1−300 M_⊙, with stars on the lower end having been ejected from the star-forming region. The large number of ~100 M_⊙ stars produced in both these simulations implies copious pair-instability supernovae should have preceded the formation of Population II stars. Susa et al. (2014) concluded that the limited resolution of their simulations likely results in a bias toward higher stellar mass, potentially by 50−100% on the upper end of the IMF. With the enormous numerical challenges of simulating the formation of the first stars, it is important to constrain theory by observations.

Despite several surveys specifically designed to discover extremely metal-poor objects in the halo (e.g., Beers et al. 1992; Cayrel et al. 2001; Christlieb 2003; Jacobson et al. 2015) and bulge (Howes et al. 2015), no metal-free stars have been detected, and only a single star indicating pollution from pair-instability supernovae (Aoki et al. 2014; but see also Aoki et al. 2015), although this may be a selection effect (Karlsson et al. 2008). A handful of extremely metal-poor stars have been identified to exhibit a chemical composition indicative of pollution from one or just a few core-collapse supernovae (e.g., Christlieb et al. 2002; Frebel et al. 2005; Caffau et al. 2011), with SMSS0313-6708 being the most likely candidate to have been polluted by a single supernova (Keller et al. 2014; Bessell et al. 2015). The remarkable nearly line-free spectrum of SMSS0313-6708indicating extremely low abundances of calcium ([Ca/H] = −7.3) and only upper limits to all heavier elements including iron ([Fe/H] < −7.5), but relatively high abundances of carbon and oxygen ([C/H] ≈ [O/H] ≈ −2.5) as well as magnesium ([Mg/H] ≈ −4.1), was found to match well the predicted yields from a massive supernova progenitor (40−60 M_⊙) and low explosion energy (1.5−1.8 B, where 1 B = 10⁵¹ erg).

Such comparisons with supernova model yields depend critically on the detailed composition, where, for example, the ratios of light to iron-peak elements and odd-even variations are the most sensitive to explosion energy and progenitor mass. Bessell et al. (2015) found that, in their analysis, the odd-even effect constrained the mass of the progenitor despite no odd elements beyond Li being detected, with the strongest constraints coming from the upper limits to N (preferring lower masses), and Na and Al (preferring higher masses). Supernova yields overall, and odd-even effects in particular, are however strongly model-dependent (see Fig. 5 of Nomoto et al. 2013), due, for example, to differences in how explosions are initiated. On the observational side, the detailed physics of spectral line formation causes deviations from LTE (this case is known as NLTE) which are known to significantly affect the Al resonance line strength. Furthermore, the interplay between effects on line formation from the steeper temperature structures of 3D model atmospheres and NLTE effects is largely unknown for all but smaller representative sets of stellar parameters (e.g., Asplund et al. 2003; Lind et al. 2013; Klevas et al. 2016), with the exceptions of Sbordone et al. (2010) and Amarsi et al. (2016a) who computed 3D NLTE grids for lithium and oxygen, respectively.

In the first analysis of SMSS0313-6708, Keller et al. (2014) used both one-dimensional (1D) hydrostatic models (ATLAS9 models from Castelli & Kurucz 2004) and ⟨ 3D ⟩temporally and horizontally averaged three-dimensional (3D) hydrodynamical models (from the Staggergrid, Magic et al. 2013a,b), with NLTE line formation for Li, Na, Mg, Ca and Fe. The analysis was improved by Bessell et al. (2015), who analyzed molecular lines using 3D LTE radiative transfer, which is appropriate and necessary due to their very strong temperature dependence (Asplund & García Pérez 2001).

The fundamental shortcoming of 1D models is that they are hydrostatic and thus must implement convective energy transfer using a parametrized mechanism; typically, mixing-length theory (Böhm-Vitense 1958). Additionally, 1D models use a free local velocity parameter called microturbulence (v_mic), representing small-scale velocity fields which serve as a line desaturation mechanism (Asplund et al. 2000b).

While ⟨ 3D ⟩temporally and horizontally averaged hydrodynamical models are more realistic than their 1D hydrostatic counterparts, the averaging process loses physical properties of the 3D atmospheric structure including surface granulation (i.e., horizontal inhomogeneities) and local velocity fields. The importance of surface granulation has been demonstrated in the past (e.g., Collet et al. 2007; Magic et al. 2013b; Dobrovolskas et al. 2013), where significant effects were found in LTE abundance analyses of metal-poor red giant stars. Magic et al. (2013b) found 3D−1D abundance corrections for weak low-excitation Fe i lines of −0.3 dex, and 3D−⟨ 3D ⟩LTE abundance corrections of −0.1 dex, indicating that the effects of surface granulation and the steeper temperature gradient have similar magnitude. In contrast, Dobrovolskas et al. (2013) found abundance corrections for weak resonance lines of Fe i to be −0.5 dex for both 3D−1D LTE and 3D−⟨ 3D ⟩LTE, indicating that the effects of granulation entirely dominate over those of the steeper temperature gradient. Using the same models and averaging techniques as Dobrovolskas et al. (2013), Klevas et al. (2016) compared 3D NLTE calculations to ⟨ 3D ⟩NLTE and 1D NLTE calculations for Li i, finding similar line strengths in 3D NLTE and 1D NLTE, in line with previous results (Asplund et al. 2003). They find that as 3D LTE line formation is inappropriate for Li i, the NLTE and 3D effects cannot simply be added together.

We present here an abundance analysis of SMSS0313-6708employing NLTE line formation in a 3D hydrodynamical model atmosphere for Li, Na, Mg, Al, Ca and Fe, using state-of-the-art atomic models with accurate collisional data. The models of Al and Fe in particular use recent accurate transition rates for hydrogen collisions, which have previously been a major source of uncertainty due to the dependence on classical models with empirical scaling factors. We introduce the spectroscopic analysis, the hydrodynamical models and the 3D NLTE modeling in Sect. 2, present the 3D NLTE effects and the resulting abundances in Sect. 3, and compare the resulting abundance determinations (Sect. 3.3) to supernova yield models and discuss the inferred properties of the progenitor and explosion in Sect. 4.

2. Methods

Abundances of Li, Na, Mg, Al, Ca and Fe were determined with state-of-the-art 3D NLTE modeling using the multi3d code (Leenaarts & Carlsson 2009), with updates as described by Amarsi et al. (2016a). We present the observational data in Sect. 2.1, introduce the 3D model atmosphere in Sect. 2.2, and give details on the 3D NLTE problem in Sect. 2.3 and our model atoms in Sect. 2.4.

2.1. Observations and spectroscopic analysis

We used existing spectra of SMSS0313-6708 taken with the MIKE spectrograph (Bernstein et al. 2003) on the 6.5 m Magellan Clay Telescope and UVES (Dekker et al. 2000) on the VLT. The UVES spectrum covers the UV region 3060−3860 Å and the optical region 4790−6800 Å at resolving power R ≈ 41 000 (see Bessell et al. 2015, for details). We extracted pipeline-reduced data from the ESO archive, and co-added the 31 exposures after applying radial velocity shifts found by cross-correlation matching. The resulting spectrum has S/N ≈ 100 in the region near 3800 Å, and S/N ≈ 200 −300 in the visual region. The resonance lines of Al i and Ca ii are not available in these spectra, and we therefore analyze these lines using the MIKE spectrum (see Keller et al. 2014) which has R ≈ 28 000 and S/N ≈ 80 in this region.

For our abundance determination of lithium, magnesium and calcium, we fit the observed line profiles. Our 3D NLTE analysis is based on interpolated synthetic 3D NLTE spectra computed in a range of abundances, and we convolve the synthetic spectra with a fitted Gaussian representing the instrumental profile and rotational broadening. For the 1D modeling, this broadening also represents large-scale (macroturbulent) velocity fields. Uncertainties are based on the statistic, where we estimate the S/N from the local continuum and typically use three free parameters (abundance, broadening and radial velocity). In the 1D abundance analysis, we used a fixed value of v_mic = 2 km s^-1. For lithium, we used the Li i 6707 Å resonance line, taking into account fine structure. For magnesium, we used the two Mg i UV lines at 3829 and 3832 Å as well as the Mg ib 5167−5183 Å triplet. We rejected the strongest UV line at 3838 Å, where the continuum placement is significantly distorted by the nearby Balmer H9 line at 3835 Å. For calcium, we used the unblended Ca ii 3933 Å line. For comparison, we also determined abundances from the 3968 Å line, which is suppressed by the nearby Balmer Hε line at 3970 Å.

For sodium, aluminium and iron, we do not detect any lines and report only upper limits. We used equivalent width uncertainties estimated using Cayrel’s formula (Cayrel 1988; Cayrel et al. 2004), , where δx is the pixel sampling. We also propagate representative uncertainties in the continuum placement into the final equivalent width uncertainty. For sodium, we use the Na iD 5895 Å line which is less affected by interstellar absorption than the neighboring Na iD 5889 Å. For aluminium, we use the unblended Al i line at 3961 Å. For iron, we stack spectra in the vicinity of the unblended Fe i lines at 3440.6, 3581.2, 3719.9, 3737.1, 3820.4 and 3859.9 Å. We estimate S/N from the nearby continuum of the stacked spectrum.

2.2. Model atmosphere

We adopt the hydrodynamical model computed by Bessell et al. (2015) using the Stagger code (Nordlund 1995; Stein & Nordlund 1998), with T_eff = 5150 K, log g = 2.2, [Fe/H] = −5. The model was computed in a cartesian geometry with 240³ grid points representing a 1145² × 458 Mm³ box in the photosphere, using tailored opacities with standard alpha enhancement ([α/Fe] = 0.4) and enhanced CNO abundances ([C/H] = [N/H] = −3.0 and [O/H] = −3.6). This composition was originally tailored to HE 0107-5240 (Christlieb et al. 2004), but is sufficiently close to the actual chemical composition of SMSS0313-6708. The simulation used 12 opacity bins selected as a function of wavelength and formation height, and treats continuum scattering using the approximation described by Collet et al. (2011). We adopt the abundances derived by Bessell et al. (2015) to compute the equation of state and background continuous opacities and generate synthetic spectra. We use the reference solar abundances from Asplund et al. (2009).

For comparison purposes, we also perform the abundance analysis using a temporally and horizontally averaged one-dimensional representation of the hydrodynamical model, denoted ⟨ 3D ⟩, as well as a 1D hydrostatic MARCS model (Gustafsson et al. 2008) computed with the same abundances by K. Eriksson. The ⟨ 3D ⟩ model used in this work has been averaged on surfaces of equal continuum optical depth τ₅₀₀, but other methods exist in the literature (see Magic et al. 2013b). We use v_mic = 2 km s^-1, and discuss this choice in Sects. 3.2 and 3.3.

Fig. 1 Temperature stratification in the hydrodynamical model shown in a gray-scale histogram as a function of the continuum optical depth. Also shown are its spatial and temporal average (pink line) and a hydrostatic (1D) MARCS model with similar stellar parameters (blue line).

The temperature distribution in the hydrodynamical model is compared to the MARCS model in Fig. 1. The outer layers (log τ₅₀₀< −2) of the hydrodynamical model are predominantly cooler than the MARCS model due to the importance of adiabatic cooling caused by gas expansion (Asplund et al. 1999). However, a significant fraction of the surface exhibits temperature inversion due to radiative reheating and mechanical compression above the intergranular lanes (see e.g. Rutten et al. 2004; Magic et al. 2013b; de la Cruz Rodríguez et al. 2015).

We estimate the effect of changes in stellar parameters on the abundance determination using MARCS models in Table 1. There, we also show the influence of using model atmospheres computed with non-standard CNO abundances, compared to the standard solar-scaled composition. As this has a small effect on the abundance determination, we judge that the difference in composition between that adopted for the model atmosphere calculation and the actual composition of the star is not likely to significantly affect results.

To lessen the computational cost we select representative snapshots from the hydrodynamical simulation sequence, and reduce these in size. We found that we could reduce the model atmospheres in size from the original 240³ grid points to retain every fourth column only and interpolate to 100 horizontal layers, giving a total of 60² × 100 grid points. In the following, we always retain 100 horizontal layers and refer to models by their number of columns, . The detailed numerical tests discussed in Sect. 3.2.2 show that these simplifications introduce negligible systematic effects, of at most 0.03 dex. The small temporal variations in strength of all lines investigated here (σ ≲ 0.02 dex) indicate that a set of just five snapshots is sufficient.

2.3. 3D NLTE calculations

NLTE radiative transfer for late-type stars involves solving the equation of statistical equilibrium to compute NLTE atomic populations, while typically neglecting feedback on the model atmosphere. Details are given in several textbooks (e.g. Rutten 2003; Hubeny & Mihalas 2014) and illustrative examples for cool stars are found in several reviews (e.g. Asplund 2005; Bergemann & Nordlander 2014). A limiting case is found in the deep photosphere, where the mean free path of photons is short compared to the gradient of local physical parameters, thus resulting in atomic populations close to LTE. In the high photosphere or chromosphere, densities are sufficiently low that non-local radiation drives higher transition rates than collisions, thus resulting in NLTE atomic populations. Spectral lines form in the intermediate layers, where radiative and collisional rates compete, and thus results are sensitive to the precise prescriptions for both.

LTE radiative transfer in 3D model atmospheres is nowadays accessible even for large grids of calculations (e.g. Allende Prieto et al. 2013; Magic et al. 2013b, 2014), but only in a few cases feasible in NLTE (e.g. Asplund et al. 2003; Sbordone et al. 2010; Lind et al. 2013; Steffen et al. 2015; Klevas et al. 2016; Amarsi et al. 2016a). For comparison, while the grid of LTE spectra computed by Allende Prieto et al. (2013) took roughly one CPU-year to compute, we exceed that when solving the statistical equilibrium of iron at a single abundance in a single 3D model snapshot. The computation time for 3D NLTE calculations scales with N_xN_yN_zN_raysN_νN_iterations, where N_x,y,z the geometrical dimensions of the atmospheric model and N_ν the number of frequencies. Here, we solve the statistical equilibrium using the 3D radiation field with N_rays = 24, distributed according to Carlson’s quadrature set A4 (Carlson 1963; Bruls et al. 1999). Additionally, several snapshots from the hydrodynamical simulation are required to sample the temporal evolution. All in all, 3D NLTE calculations turn out to be some 10⁵ times more expensive than in 1D for a given model atom. In order to make these calculations feasible, care must thus be taken to reduce the complexity of both model atoms and atmospheres.

Implementation details on the version of multi3d used here are given by Leenaarts & Carlsson (2009) and Amarsi et al. (2016a). We note that the code implements isotropic coherent scattering in background processes, considering Rayleigh scattering by H i in the far red wing of Lyman-α and by H₂, as well as Thomson scattering on electrons. Background line opacities, mainly due to hydrogen, were included as true absorption.

2.4. Model atoms

Fig. 2 Term diagrams of the atoms analyzed in this work, illustrating all bound-bound transitions considered in the solution of the statistical equilibrium. Thicker lines are due to overlapping transitions. Lines or multiplets used in the abundance analysis are shown in red and labeled. The excited states of Mg ii, Al ii and Fe ii are omitted from the figure.

Due to the high computational cost of 3D NLTE calculations, the very detailed general-purpose atoms with (tens of) thousands of spectral lines typically used in 1D analyses must be modified before use. In order to lessen the computational cost, we reduce the number of radiative transitions by discarding those which do not have signficant impact on the statistical equilibrium, using an approach described by Lind et al. (in prep.). In summary, we quantify the relative importance of radiative transitions by comparing their radiative bracket, | n_ir_ij−n_jr_ji |, where r_ij are radiative rate coefficients, and removing those transitions where this value (essentially the influence on the statistical equilibrium) is relatively small. We perform this test at τ₅₀₀ = 0.1 using a representative ⟨ 3D ⟩ Stagger model, requiring that NLTE populations of states involved in the lines we investigate here remain unaffected to within 2% (0.01 dex) at depths log τ > −2, or farther out in the case of stronger lines forming in the UV. It is also required that the inferred abundance from all lines used in the analysis remains unaffected to within 0.01 dex. We apply the same criterion when determining which atomic levels to retain, and always compare the final atom to the original one. The complexity of our resulting model atoms is summarized in Table 2 and the corresponding term diagrams are shown in Fig. 2. In particular, we note that we were able to discard all bound-free transitions of the ionized species, as these did not significantly affect the populations of the relevant atomic levels, even in Ca ii. We could also significantly reduce the number of bound-free transitions in the neutral species of atoms where just a few of these were dominant.

Collisional transition rates for impact by electrons and hydrogen atoms are a major source of uncertainty for NLTE modeling. Neutral hydrogen atoms are orders of magnitude more numerous than free electrons, but the low mass of electrons leads to high thermal velocities and thus large impact rates on atoms. Hydrogen atoms have low impact rates as they are relatively massive and thus slow, and according to the Massey criterion (Massey 1949) this also leads to the expectation of adiabatic collisions with smaller inelastic collisional cross-sections. As calculations of accurate hydrogen collisional rates have become feasible, it is apparent that both hydrogen and electron collisions are important, and that classical recipes describing hydrogen collisions (Drawin 1968) are typically far worse than order-of-magnitude estimates (see Barklem et al. 2011; Barklem 2016b, for recent reviews on the subject). These inaccuracies have traditionally been compensated for, to some degree, by introducing empirical scaling factors, which vary significantly from element to element as well as between authors (see e.g.Table 2 of Bergemann & Nordlander 2014). The model atoms used in this work use quantum mechanical calculations for hydrogen transition rates, and are free from empirical scaling factors.

2.4.1. Lithium and sodium

We use the same atoms as Lind et al. (2013), meaning slightly modified versions of those presented by Lind et al. (2009a), Lind et al. (2011). The Na atom is truncated at n = 10, leaving 20 levels of Na i, and for both atoms we consider only the 91 most important radiative transitions. Hydrogen collision transition rates (Barklem et al. 2003, 2010) are based on quantum scattering calculations based on quantum chemical data for the LiH and NaH molecules (Belyaev & Barklem 2003; Belyaev et al. 2010). For lithium, the dominant electron collisional rates come from Park (1971). While the R-matrix calculations performed by Osorio et al. (2011) differed significantly, the dominant transitions happen to agree sufficiently well that the effects on the resonance line strength are smaller than 0.01 dex.

2.4.2. Magnesium

We use a simplified but realistic Mg model atom assembled for this work. We adopt hydrogen collision rates from Barklem et al. (2012), based on quantum mechanical cross-sections (Guitou et al. 2011; Belyaev et al. 2012) for low-lying states of Mg i, based on quantum chemical data for the MgH molecule from Guitou et al. (2010). For Rydberg transitions we computed rates based on the free electron model of Kaulakys (1985, 1991) using the implementation of Barklem (2016a). The relative importance of electron and hydrogen collisional transitions was demonstrated by Osorio et al. (2015), indicating the strong sensitivity to hydrogen collisional rates in giants. We adopt energy levels, f-values, and photoionization cross-sections from TOPbase (Cunto et al. 1993; Butler et al. 1993; Mendoza et al. 1995), except for the Mg ib lines where we adopt the laboratory f-values from Aldenius et al. (2007). The TOPbase photoionization cross-sections are tabulated using a dense frequency sampling necessary to resolve narrow resonance features. However, limitations in the calculations result in minor uncertainties in the positions (but not necessarily strengths) of resonances, which may spuriously coincide with narrow background opacity features. To reduce the potential bias of such coincidences, Bautista et al. (1998) suggested smoothing photoionization cross-sections to represent the inherent uncertainty of resonance energies. Allende Prieto et al. (2003) executed such smoothing, and we use their data, interpolated onto an equidistant frequency grid in order to sample the frequency dependence of background opacities.

Our Mg atom is represented by 43 levels of Mg i, 2 levels of Mg ii and the Mg iii continuum. For Mg i, we retain fine structure in the 3p ³P^o term, and l-splitting when l ≤ 3 and n ≤ 6. Levels with a given principal quantum number n and l ≥ 4 or n> 6 were combined, regardless of spin. In line with Osorio et al. (2015), we found that the collisional couplings were sensitive to the inclusion of high-lying states, and thus include states all the way up to n = 80 in a coarse set of superlevels representing states with n = 9−10, 11−12, 13−14, 15−20, 21−30, 31−50, 51−80. For Mg i, we include 70 spectral lines, 47 of which originate from the 3p ³P^o states, and continua coupling the Mg ii ground state to the 3p ³P^o and 3p ¹P^o states with photoionization edges at 2512 and 3756 Å (driving over-ionization), and five higher excited states with edges at 2.2 μm and beyond (driving over-recombination). The two levels of Mg ii, 3s ²S and 3p ²P^o, are collisionally and radiatively coupled via the resonance line at 2799 Å. In total, the bound-bound and bound-free radiative transitions are computed at 1600 frequency points.

2.4.3. Aluminium

We use a model atom to be described by Nordlander et al. (in prep.), with hydrogen collisional rates from Belyaev (2013a) computed for low-lying states using a model approach based on semi-empirical couplings and the Landau-Zener model (Belyaev 2013b), and rates computed following Kaulakys (1985, 1991) for Rydberg states. Electron collisional rates come from open-ADAS¹ (computed as described by Badnell 2011) when possible, including forbidden transitions, and are otherwise computed using Seaton (1962) for allowed transitions. Energy levels and f-values are taken from NIST (Kelleher & Podobedova 2008) when possible, or otherwise the most recent data from Kurucz (Kurucz & Bell 1995). Photoionization cross-sections are again taken from Allende Prieto et al. (2003), based on data from TOPbase (Mendoza et al. 1995).

Our Al atom is represented by 42 levels of Al i and 2 levels of Al ii. For Al i, we resolve fine structure in the 3p ²P^o ground state, and we retain l-splitting when n ≤ 8, and otherwise represent the states by a single superlevel at a given n-value when 9 ≤ n ≤ 15. States with 16 ≤ n ≤ 20 are represented by a single superlevel. Levels of Al i are coupled radiatively by 90 b-b transitions, 40 of which are resonance lines. We also include four bound-free transitions, coupling the fine-structure 3p ²P^o ground states and the 3d ²D and 4p ²P^o states to the Al ii continuum with ionization edges at 2075, 6312 and 6528 Å, driving over-ionization. The bound-bound and bound-free transitions use a total of 2100 frequency points. Al ii is represented by the ground and first excited state, with only collisional coupling.

2.4.4. Calcium

We modified the model atom of Lind et al. (2013). Only lines of Ca ii are visible in the spectrum of SMSS0313-6708, and as the dominant ionization stage, the behavior of Ca ii does not noticeably depend on Ca i. We thus only retain the 4s ²S ground state of Ca ii plus the first two excited states (3d ²D and 4p ²P^o, with resolved fine-structure) and the Ca iii continuum. Electron collisional rates are based on R-matrix calculations (Meléndez et al. 2007), but the statistical equilibrium is not very sensitive to collisional rates. We use energy levels from TOPbase (HE Saraph & PJ Storey, to be published) and f-values from Kurucz (Kurucz & Bell 1995) which are very similar to the data in NIST.

By running tests on a complete atom representing Ca i by 39 levels, 244 b-b transitions and 39 b-f transitions, we verified that NLTE effects on Ca i indeed do not significantly affect the lines of Ca ii on the level of 0.01 dex. For these tests, we included hydrogen collisions for Ca i based on Drawin (1968) rescaled by S_H = 0.1 (following Mashonkina et al. 2007), as Lind et al. (2013) found acceptable fulfillment of the excitation and ionization equilibria using these collisional rates for a set of very metal-poor dwarfs and subgiants. We note that Barklem (2016c) recently computed rates for charge transfer and collisional excitation of Ca i, but do not expect these transition rates to affect the conclusion that the statistical equilibrium of Ca ii is not significantly affected by NLTE effects on Ca i.

We model the resonance Ca ii H and K lines at 3933−3968 Å and the Ca ii IR triplet at 8500 Å, using 550 frequency points. We neglect continuous transitions coupling Ca ii to the Ca iii continuum as the photoionization cross sections are relatively small, with edges at λ < 1420 Å, resulting in negligible influence on the statistical equilibrium.

2.4.5. Iron

We utilize a heavily reduced version of the newly compiled model atom which will be described in an upcoming paper (Lind et al., in prep.), and has previously been applied by Amarsi et al. (2016b). This model atom adopts collision rates for Fe+H collisions calculated with the asymptotic two-electron method presented by Barklem (2016c) and applied there to Ca+H. The calculation used here includes 138 states of Fe i, and 11 cores of Fe ii, leading to the consideration of 17 symmetries of the FeH molecule. This removes the systematic uncertainty of classical collisional rates (Drawin 1968), which in previous works were rescaled by a factor S_H calibrated on standard stars. Electron collisional rates were computed using the formula from van Regemorter (1962). The statistical equilibrium of iron represents a formidable problem, as the atom is vastly more complex than the others investigated in this work. Even disregarding fine structure splitting, the full model atom considers 10⁵ radiative transitions, requiring sampling at 10⁶ frequency points, at an estimated computational cost of several millennia of CPU time for the full 3D NLTE solution.

Our reduced atom is represented by 421 levels of Fe i, 12 levels of Fe ii, and the Fe iii continuum. We have reduced the number of Fe ii levels to retain only the nine low-lying states required to correctly compute the partition function, as well as the upper states of a few resonance multiplets in the UV (z ⁶P^o, z ⁶D^o and z ⁶F^o). Additionally, we do not resolve fine structure when solving the statistical equilibrium, and include only 776 and 12 spectral lines for Fe i and Fe ii, respectively, as well as 33 continuous transitions coupling Fe i to Fe ii. In total, these radiative transitions are computed at 12 000 frequencies.

3. Results

We discuss the NLTE mechanisms and the resulting 3D NLTE effects in Sect. 3.1 and in Sect. 3.2 where we also compare 1D modeling. We estimate errors in Sect. 3.2.2, and present the results of the abundance analysis in Sect. 3.3.

3.1. Resolved disk-center 3D NLTE effects

Fig. 3 Spatially resolved NLTE equivalent width correction at disk-center intensity (μ = 1), shown for representative lines of Li i (6707 Å, top left), Al i (3961 Å, top right), Ca ii (3933 Å, bottom left), and Fe i (3719 Å, bottom right) in a high-resolution snapshot of the hydrodynamical model (2402 grid points for Li, Al and Ca; 1202 for Fe). Although the behavior is shown for vertical rays, we emphasize that the statistical equilibrium was solved using the 3D radiation field. The NLTE line strength correction, Δlog Wλ,μ = 1 = log Wλ,μ = 1(NLTE)−log Wλ,μ = 1(LTE), corresponds roughly to the abundance correction but with opposite sign. Red and blue colors indicate weaker and stronger lines under NLTE than under LTE, meaning positive and negative abundance corrections. For this illustration, we adopt the abundances A(Li) = 0.70, A(Al) = 1.15, A(Ca) = −0.62, and A(Fe) = 0.60. Na, not shown, behaves very similarly to Li, and the corresponding illustration for Mg is shown in the bottom left of Fig. 4.

We illustrate the predicted NLTE line strength corrections at disk-center (i.e., using vertical rays, where μ = 1) on the resolved surface of a high-resolution snapshot from the hydrodynamical simulation for representative lines of Li i, Al i, Ca ii and Fe i in Fig. 3. The NLTE line strength correction, Δlog W_{λ,μ = 1} = log W_{λ,μ = 1}(NLTE)−log W_{λ,μ = 1}(LTE), represents the difference in line strength between the NLTE and the LTE case, and is approximately the same as the abundance correction (but with opposite sign) for weak lines. Note, however, that line strengths vary from center to limb with typically stronger NLTE−LTE differences toward the limb as rays with small μ probe higher atmospheric layers where collisional rates are lower.

Below, we give detailed descriptions of the NLTE behavior of each species at disk-center. We find the NLTE corrections across the surface of the simulation to be negative across the bulk of the hot granules. The overall behavior varies between elements due to different processes being responsible for NLTE effects, with our different diagnostic lines forming from levels of different excitation energy as well as ionization stages.

3.1.1. Lithium and sodium

At low abundances, the Li and Na resonance lines behave similarly as they form in the transition region between over-ionization and over-recombination. Over-ionization tends to dominate in granules while over-recombination gives a net line strengthening in the intergranular lanes. The interplay depends on the line strength, as the ground state tends to be over-populated in deeper layers but depleted in higher layers. Since stronger lines form in higher layers, the NLTE correction becomes more negative the higher the abundance. As long as the lines are very weak, this ionization balance effect is more important than photon losses in the core.

3.1.2. Magnesium

Fig. 4 Spatially resolved line equivalent widths of the Mg i 5183 Å  line at disk-center intensity (μ = 1) in a high-resolution (2402) snapshot of the hydrodynamical model, assuming A(Mg) = 3.80. Top left: NLTE (right: LTE) equivalent widths. Bottom left: NLTE equivalent width correction. Red and blue colors carry the same meaning as in Fig. 3. Bottom right: granulation pattern imprinted on the relative continuum intensity.

We illustrate the predicted line strength for the Mg ib 5183 Å line on the resolved surface of the hydrodynamical simulation in Fig. 4. In LTE, the predicted line strength correlates with the continuum intensity I_c, that is, the optical granulation pattern, such that the line is typically stronger in the bright, granular upflows than in the cooler intergranular lanes. This is because the saturated line core forms high above the continuum, where temperature inversions above intergranular lanes reduce line strengths. The NLTE line strength in the granules is however clearly suppressed due to over-excitation driven mainly by the UV triplet at 3829−3838 Å, which efficiently depopulates the 3p ³P^o states. The fact that the NLTE effect is strongest in the hot regions is expected, as over-excitation and over-ionization result from the J_ν−B_ν split, which is larger where the temperature structure is steeper. Consequently, this effect is typically smaller in the intergranular lanes, with cooler regions even experiencing enhanced line strength in patches with significant temperature inversions.

3.1.3. Aluminium

The ground state of Al is efficiently de-populated relative to LTE by over-ionization and over-excitation driven by the resonance lines. The NLTE line strength correction varies strongly between essentially zero in the cool intergranular lanes and as much as −0.7 dex in the granules, where the strong temperature gradient results in a large J_ν−B_ν, which drives the NLTE effect.

3.1.4. Calcium

For Ca, the dominant NLTE driver is over-excitation of the upper state, 4p ²P^o, driven by the resonance lines. This has minute effects on the ground state but strongly over-populates the upper states. A flow back down via the IR triplet lines over-populates the 3d ²D state as well, which is collisionally- but not radiatively coupled to the ground state. The resonance line formation behaves like a two-level atom, with a superthermal source function (due to the population inversion) equal to the mean intensity , resulting in weakened lines. Aside from the same type of saturation effect as was seen for Mg in Fig. 4, in regions where the temperature structure is inverted in layers high above the continuum, NLTE effects on Ca ii vary only weakly across the surface.

3.1.5. Iron

NLTE effects on the Fe lines are similar to those seen for Al. The same drivers; photon pumping in the resonance lines and over-ionization, are responsible for the departures from LTE. However, while photoionization operates on the ground state of Al i, it is mainly the excited states of Fe i which are over-ionized. NLTE effects on the Fe i resonance lines are thus due to the resonance lines driving a flow from the ground state to the intermediate states, which are over-ionized. The NLTE line strength correction is stronger (more negative) than for Al, with a broad distribution of line strength corrections peaking at −1 dex but ranging to −1.5 dex.

3.2. 3D NLTE effects on disk-integrated flux spectra

Line strengths computed in NLTE and LTE using the highest resolution (240², except for Fe where 120² was used) snapshot from the hydrodynamical simulation are indicated in Fig. 5 by horizontal dashed lines. NLTE line strengths are smaller (giving positive abundance corrections), with differences ranging from as little as 0.04 dex for the Na i 5895 Å line to 1.15 dex for Fe i 3719 Å. We note, and caution, that the NLTE effects vary between lines of the same species, as well as with abundance, and cannot necessarily be directly applied to an abundance analysis. Additionally, the lines of Mg and Ca are saturated, leading to larger effects on abundances than what is indicated by differences in line strengths.

In the following subsections, we compare these 3D NLTE results to those computed using 1D or ⟨ 3D ⟩ models and to calculations using model atmospheres reduced in size, as well as other snapshots from the hydrodynamical simulation, and evaluate the influence of non-vertical radiative transfer on the NLTE solution.

Fig. 5 Variation of predicted disk-integrated line strengths, as inferred from a snapshot from the hydrodynamical simulation downsampled from the original 2402 grid points as well as the horizontally averaged ⟨ 3D ⟩ and hydrostatic 1D case. A vertical dashed line separates the 3D results from the 1D hydrostatic and ⟨ 3D ⟩ results. Calculations in NLTE (black diamonds) and LTE (blue squares) are illustrated for representative lines using A(Li) = 0.70, A(Na) = 0.89, A(Mg) = 3.80, A(Al) = 1.15, A(Ca) = −0.62, and A(Fe) = 0.60. When additional snapshots have been computed (at resolution 602), their individual results are indicated by horizontal lines. A red greek cross indicates the corresponding NLTE result computed in the so-called 1.5D approximation (at resolution 602). Each panel spans over 0.35 dex, except for Al and Fe due to their very strong NLTE effects. For clarity, the line strength determined from the highest-resolution model is indicated as a horizontal dashed line across each panel, and representative errors of 0.03, 0.06 and 0.09 dex are indicated by gray shaded regions.

3.2.1. 3D−1D comparison

Line strengths computed using ⟨ 3D ⟩ horizontally and temporally averaged one-dimensional structure, and the 1D hydrostatic MARCS model, are shown in Fig. 5. Comparing LTE results, we find that ⟨ 3D ⟩ models reproduce the 3D LTE results well, with deviations of less than 0.05 dex except for Fe where lines are too weak by 0.2 dex. In ⟨ 3D ⟩ NLTE, lines of Li, Na, Al and Fe are too strong by 0.04−0.07 dex compared to 3D NLTE. The Ca line is too strong by 0.11 dex, and the lines of Mg are too strong by 0.12−0.17 dex, depending on the line.

For MARCS models in LTE, differences with respect to 3D LTE are larger than those for ⟨ 3D ⟩ models, as expected, with Fe line strengths being 0.3 dex too weak. In NLTE, results for MARCS models are very close to 3D NLTE for Li, Na and Al, to within 0.02 dex. Lines are too strong by 0.05 dex for Fe, 0.10 dex for Ca, and 0.10−0.16 dex for Mg.

Overall, NLTE predictions for weak lines using these 1D MARCS or ⟨ 3D ⟩ models are not too different from the full 3D NLTE solution. For the saturated lines of Mg and Ca however, line strengths are significantly over-estimated, indicating that the adopted value v_mic = 2 km s^-1, may be too high. As the curve of growth is flat for saturated lines, the resulting effect on abundances is even larger than that indicated by the difference in line strength. The choice of v_mic parameter is discussed further in Sect. 3.3.

3.2.2. 3D modeling errors

We investigated both the influence of reducing the resolution of a given snapshot from the hydrodynamical simulation, and the temporal snapshot-to-snapshot variations, from data illustrated in Fig. 5.

Reducing the resolution appears to introduce systematic errors in both NLTE and LTE calculations. NLTE results for lower-resolution models systematically overestimate line strengths for Li, Na, Mg and Ca, while underestimating them for Al and Fe. The magnitude of these effects appears larger for lines where NLTE effects are large. Errors are smaller in LTE, and are similar for all lines.

We find that for 3D NLTE modeling, the model atmospheres can safely be reduced in resolution to 60², with errors less than 0.03 dex compared to calculations at the maximum resolution. For Fe i, we have not been able to solve the statistical equilibrium for models at resolutions higher than 120². Line strengths deduced from these models differ from results at resolution 60² by only 0.02 dex. As the NLTE effects are similar to those seen for Al i, albeit stronger, we estimate that the errors introduced by running these simulations at lower resolution should not be larger than ~0.03 dex.

Little work has been carried out to test the sensitivity of abundance analyses to the numerical resolution used in the 3D stellar models and in the radiative transfer calculations. Asplund et al. (2000a) studied iron lines in the solar spectrum, where the low-resolution models (50² compared to 200²) were 50−100 K warmer in the outer atmospheric layers. In their LTE synthesis, this resulted in weaker lines indicating higher abundances by ~0.05 dex, but reaching as high as 0.1 dex for low-excitation lines. Systematic temperature differences of similar magnitude have been found at higher resolution (140² × 150 compared to 400² × 300) in simulations of the solar photosphere computed with the CO⁵BOLD code (Freytag et al. 2012). Based on the small effects seen at very high resolution (480²) in models computed for HD 122563 by Collet et al. (2009), we expect that our hydrodynamical model is sufficiently large in this respect.

In addition to the numerical resolution, we analyzed sequences of snapshots to estimate the effect of temporal variations, and indicate the line strengths inferred from individual snapshots by horizontal lines in Fig. 5. The snapshot-to-snapshot RMS variation at 60² resolution is less than 0.02 dex. The set of five snapshots we use for the abundance analysis does sample the simulation well, and even in the case of Fe i, where we use a single snapshot, the small sample of snapshots does not introduce significant errors.

In summary, we find that it is unlikely that combined errors due to the limited resolution and number of snapshots is larger than 0.05 dex for any species analyzed here, as compared to the full 240² model. We find smaller errors for Li i, Na i, Mg i and Ca ii, with combined effects likely not larger than 0.02 dex. The 3D and NLTE effects are thus much greater than these numerical errors. Additionally, these errors fall well within the estimated error bars associated with uncertainties in the stellar parameters, indicated in Table 1.

Other potential sources of errors include, for example, effects of the opacity binning and NLTE effects on hydrogen which may affect the atmospheric structure, and input data such as photoionization and collisional transition rates in the NLTE modeling. Rigorous propagations of atomic data uncertainties are beyond the scope of this analysis, but see Osorio et al. (2015), for example, where the influence of different collisional processes was explicitly tested for Mg.

3.2.3. The influence of non-vertical radiative transfer

To help gauge the influence of non-vertical radiative transfer on the NLTE solution, we also performed test calculations using the so-called 1.5D approximation, where the NLTE solution is acquired separately, column-by-column in the model atmosphere. The 1.5D approximation can also be utilized to save computational resources, as a full 3D solution of the statistical equilibrium is avoided (e.g., Kiselman 1993; Shchukina & Trujillo Bueno 2001; Shchukina et al. 2005, 2009; Pereira & Uitenbroek 2015). In the case of the 777 nm oxygen triplet lines, Amarsi et al. (2016a) found that solving the statistical equilibrium using the full 3D radiation field did not differ greatly from the 1.5D approximation, but this is not a general result applicable to any atom or type of star.

Line strengths computed in the 1.5D approximation are indicated in Fig. 5. We find that two competing mechanisms appear to be important. With the full 3D radiation field, regions above the cool intergranular lanes are illuminated by the nearby hot granules via inclined rays. As these rays are missing in 1.5D, the J_ν−B_ν split decreases, thus leading to strengthened lines in these regions. Conversely, the radiation field above the hot granules is somewhat stronger in 1.5D, leading to weakened lines in these regions. We find that the net result is a slightly more positive NLTE line strength correction (leading to lower abundances) for all six elements, albeit with effects of at most 0.02 dex for Li, Na and Ca, 0.04 dex for Mg, 0.10 dex for Al, and 0.22 dex for Fe.

3.3. Abundance analysis

Fig. 6 Upper limit determination for iron, comparing stacked spectra centered on the six strongest unblended lines, Fe i 3440.6, 3581.2, 3719.9, 3737.1, 3820.4 and 3859.9 Å. The 3D NLTE line profiles are shown for the 1σ, 2σ and 3σ upper limits, corresponding to line strengths of 0.7, 1.5 and 2.2 mÅ. The small bump at + 0.4 Å in the synthetic spectrum is due to the nearby line Fe i 3441.0 Å.

The results of our 3D NLTE analysis are shown in Table 3 and illustrated for iron in Fig. 6. The stacked spectrum of the Fe i 3440.6, 3581.2, 3719.9, 3737.1, 3820.4 and 3859.9 Å lines does hint at a weak absorption feature, but we stress that it is consistent with the noise level at 1σ. For comparison, we also report in Table 3 the corresponding LTE results, as well as results from one-dimensional modeling using ⟨ 3D ⟩and 1D MARCS models. We combine our results with abundances of C, N and O based on 3D LTE molecular line analyses, and 1D LTE analyses for the remaining elements from Bessell et al. (2015) in Table 4. These are our recommended abundances for comparisons to supernova models.

Differences between our 3D NLTE results and abundances determined in LTE using 1D MARCS models are negligible for Li and Na, but amount to + 0.5 dex for Mg, Al and Ca, and + 0.8 dex for Fe. Comparing instead to 1D NLTE results, we find excellent agreement for Li, Na and Al, and differences decrease somewhat to + 0.2 dex for Mg and Fe, and + 0.3 dex for Ca. In contrast, abundances based on the 3D LTE analysis deviate more strongly from 3D NLTE than did the 1D LTE analysis, indicating that 3D LTE analyses of atomic lines are unsuitable in the presence of strong NLTE effects. This is expected, as the steeper temperature structures of the 3D models give rise to even stronger NLTE effects (as is evident from Table 3 for every element except Na).

Abundance analyses in LTE using the temporally and horizontally averaged ⟨ 3D ⟩ model, which has a steeper temperature gradient than the 1D MARCS models, again deviate more strongly from the 3D NLTE results than did the 1D LTE analysis, and are more similar to 3D LTE. In NLTE, results with the ⟨ 3D ⟩ model are similar to the 1D NLTE analysis, indicating that adopting the more realistic temperature structure of the ⟨ 3D ⟩ model is not sufficient, but in fact horizontal inhomogeneities, granulation, must be taken into account to correctly model the NLTE line formation. Additionally, results for ⟨ 3D ⟩ models are sensitive to the averaging method (Magic et al. 2013b).

Fig. 7 Abundance results for five lines of magnesium using different types of modeling, shown as a function of the reduced equivalent width. Line central wavelengths are indicated in units of Å. For the 3D NLTE and 3D LTE cases, we overplot linear fits to the trends. In 3D NLTE, we find a positive slope of less than 2σ significance, while 3D LTE as well as the ⟨ 3D ⟩ results indicate >3σ significant negative slopes.

As illustrated in Fig. 7, Mg abundances derived from the individual lines in 3D NLTE are in excellent agreement, with a scatter of only 0.03 dex, weakly correlated with line strength. 3D LTE results, on the other hand, indicate significantly lower abundances derived from the stronger lines. ⟨ 3D ⟩ NLTE and LTE abundances also strongly anticorrelate with line strength; we find very similar behavior using 1D MARCS models. As all lines but the weakest are saturated, 1D and ⟨ 3D ⟩results are sensitive to the choice of the microturbulence parameter, which we set to v_mic = 2 km s^-1. Decreasing v_mic by 0.5 km s^-1 would increase the abundance determined from the strongest line by 0.10 dex. Thus, in order to bring lines of Mg i of different strength into agreement, v_mic would need to decrease significantly, to 1 km s^-1. Such a low value would increase the average deduced Mg abundance by 0.1 dex, bringing ⟨ 3D ⟩ NLTE results into agreement with 3D NLTE. Setting v_mic = 1 km s^-1 would however also weaken the heavily saturated Ca ii 3933 Å line by 0.1 dex, causing an enormous 0.6 dex increase in the inferred abundance. While we do not use the Ca ii 3968 Å line in our abundance analysis, we note that results from the two lines agree to within 0.01 dex in 3D NLTE. Setting v_mic = 1 km s^-1 would lead to a 0.2 dex disagreement between the two Ca ii lines in ⟨ 3D ⟩ NLTE. Therefore, the free parameter v_mic cannot be invoked to simultaneously flatten trends with line strength of both Mg i and Ca ii, nor to bring 3D NLTE and ⟨ 3D ⟩ NLTE synthesis of these saturated lines into agreement.

Compared to the results of Keller et al. (2014) and Bessell et al. (2015), our abundances are typically slightly higher when running comparable analysis methods. To explain this, we first note that our tailored atmospheric models were computed using a chemical composition similar to SMSS0313-6708, affecting the overall structure with somewhat higher T_eff (by 25 K) and lower log g (by 0.1 dex) than was determined by Keller et al. (2014). Our abundance analysis is based on spectra computed with multi3dusing the hydrodynamical model. In contrast, Keller et al. (2014) and Bessell et al. (2015) used 1D models interpolated to the adopted stellar parameters in their 1D LTE analyses and then applied representative abundance corrections based on comparisons between 1D LTE and 3D LTE or ⟨ 3D ⟩ NLTE modeling. While differences in stellar parameters are well within the estimated observational error bars of 100 K and 0.2 dex, they lead to us determining higher abundances by 0.02−0.04 dex, or lower abundances by 0.05 dex for Ca, based on the sensitivities listed in Table 1. Second, the UVES spectra were aligned and co-added anew. Upper limits to the faint features of Al, Na and Fe are determined using Cayrel’s formula (Cayrel 1988; Cayrel et al. 2004) taking into account uncertainties in the continuum placement, and for Fe we stack a larger number of spectral lines.

Our ⟨ 3D ⟩NLTE abundances of Li, Na, Mg and Al agree to within 0.1 dex with Bessell et al. (2015). Discrepancies are likely due to differences in the data reduction, continuum normalization, and, in the case of Na and Al, the upper limit determinations. For Mg, our 3D NLTE abundance is higher than the previous best estimate by 0.25 dex. For Al, our 3D NLTE abundance limit (1σ) is higher than the previous 1D LTE value of Bessell et al. (2015) by 0.5 dex. This is entirely due to the 3D NLTE effect, as our 1D LTE 1σ upper limits agree to within 0.02 dex. Our ⟨ 3D ⟩NLTE Ca abundance agrees with previous estimates, and our 3D NLTE abundance is higher than this by 0.3 dex.

Finally, our 3σ upper limit for Fe is A(Fe) = 0.97 in 3D NLTE, compared to A(Fe) = −0.02 as determined in ⟨ 3D ⟩NLTE by Bessell et al. (2015). Our higher abundance limit is mainly due to our more conservative statistical approach to estimating the upper limit (+ 0.5 dex), new NLTE data (+ 0.2 dex) and the use of 3D NLTE rather than ⟨ 3D ⟩NLTE (+ 0.1 dex), all of which lead to a higher abundance. This does not however fully explain our higher abundance limit, indicating that additional effects due to differences in, for example, synthesis codes, averaging methodology and continuous opacities may be present. Tracing and quantifying these differences is challenging, but we note that our analysis is directly based on 3D NLTE synthesis, rather than 1D LTE synthesis with a posteriori corrections. We stress that at A(Fe) = 0.97, the stacked 3D NLTE spectrum has an equivalent width of just 2.2 mÅ, and that reducing the abundance by 0.5 dex would result in a line strength of 0.7 mÅ.

4. Discussion

Fig. 8 Best fitting core-collapse supernova model yields (Heger & Woosley Heger & Woosley 2010) mixed with primordial BBN-composition gas (see text), compared to the inferred abundances and upper limits listed in Table 4. Labeled elements are based on 3D NLTE analysis of atomic lines (blue characters, this work) or 3D LTE analysis of molecular lines (orange characters, from Bessell et al. 2015). Other upper limits are based on 1D LTE analyses (from Bessell et al. 2015). The progenitor stellar mass, explosion energy, and amount of mixing of the best-fit model are indicated in the legend.

We compare our recommended abundances (Table 4) to core-collapse supernova yields computed for Population III stars (Heger & Woosley 2010, and subsequent online updates in 2012²) in Fig. 8. The recommended abundances are based on our 3D NLTE synthesis for Li, Na, Mg, Al, Ca and Fe. We use literature data from Bessell et al. (2015), where the abundances of C, N and O are based on 3D LTE synthesis of molecular features. Other elements, including most heavier elements (Z> 20), have only upper limits derived from 1D LTE analyses. As these upper limits are significantly higher than the supernova yields, we do not expect errors due to unknown 3D NLTE effects on these elements to significantly bias our results.

Fig. 9 Best fitting supernova model yields assuming a selection of different progenitor masses, where thicker lines and symbols represent higher mass. The explosion energy (varying color) and mixing fraction (varying plot symbol) have been selected as the best fit to the observed abundances at each selected mass. Residuals are shown beneath for all available elements, with the 1σ, 2σ and 3σ measurement uncertainties indicated in dark to light shades of gray.

We consider a scenario where the supernova ejecta enrich primordial gas containing He and Li produced in Big Bang Nucleosynthesis (abundances from Cyburt et al. 2016, i.e., A(Li) = 2.67 ± 0.06 and Y_p = 0.247). The final lithium abundance is then reduced according to the typical difference between observed abundances of Li in metal-poor low-mass stars on the lower RGB, A(Li) = 1.1, (e.g., Lind et al. 2009b) and the BBN level, a difference of 1.6 dex. As it is not possible to accurately predict the amount of intrinsic Li depletion in extremely metal-poor stars (see, e.g., Sbordone et al. 2010; Meléndez et al. 2010), we do not use Li in the model constraint. The shown yield for H indicates the logarithmic enrichment factor, that is, the logarithmic ratio of enriched to primordial hydrogen. This represents the dilution of the supernova ejecta, which is a free parameter set to reproduce the observed abundance levels.

We compare yields to observations using the starfitcode (Heger & Woosley 2010). The yields were computed for stellar evolution models of non-rotating metal-free stars spanning a wide range in mass (9.6−100 M_⊙). As a robust physical explosion mechanism is not yet known, the models are “exploded” by means of a piston action, exploring a range of discrete energies (0.3−10 B, where 1 B = 10⁵¹ erg). Additionally, these one-dimensional models cannot predict the amount of Rayleigh-Taylor induced mixing, and the models are instead mixed by means of an averaging process parametrized as a fraction of mass relative to the He core, f_mix. As yields vary smoothly as a function of energy, albeit sometimes rapidly, we interpolate the grids to a more densely spaced set of explosion energies before running starfit. Variations as a function of mass are however rather chaotic and cannot be considered continuous (clearly illustrated in Fig. 10 of Heger & Woosley 2010).

The best-fitting supernova model has a progenitor of intermediate mass (37 M_⊙), exploding with essentially no mixing (log f_mix ≤ −2.2, or 0.6%) and low energy (2 B). The low amount of mixing is consistent with predictions from hydrodynamical simulations (Joggerst et al. 2009). Increasing the amount of mixing leads to higher yields of Na, Al and Si, and models with log f_mix> −2.2 produce yields in disagreement with our upper limits. Due to the low explosion energy, material near the core is not imparted sufficient kinetic energy to leave the gravitational well, but instead falls back into the nascent black hole. The explosion energy thus correlates with yields of heavy elements including O but not C. Increasing the explosion energy by just 10% produces [O/C] = 0.6 and [Ca/H] > −6, in strong disagreement with observations.

The sensitivity to mass is less straightforward, as we find that a wide range of models can reproduce observations well, but only if the explosion energy is allowed to vary. We illustrate the yields from a selection of well-fitting models in Fig. 9. The models have been selected as those which best match observations for any explosion energy and mixing efficiency within a range of progenitor masses. Hypernovae with explosion energies at 5 B do not match observations as they over-produce Na and Al while underproducing C and O compared to Mg. At higher energies, hypernovae also strongly overproduce calcium and iron-peak elements.

More massive stars in the initial mass range 140−260 M_⊙ are expected to explode as pair-instability supernovae, with a well-understood explosion mechanism. In the grid by Heger & Woosley (2002), models show a strong odd-even effect (typical of such supernovae) similar to our observations, but strongly overproduce heavier elements compared to C and O.

Fig. 10 Comparison of progenitor mass and explosion energy among the well fitted (χ2< 3) supernova models. Colors (explosion energy) and symbols (mixing) carry the same meaning as in Fig. 9. Results are shown on the left for our recommended abundances (from Table 4) and on the right using 1D LTE abundances. To guide the eye, power-laws fitted to the low-mass and high-mass models in the left panel are shown also to the right.

The starfitcode selects models using a χ² statistic calculated according to Eq. (4) of Heger & Woosley (2010), where we constrain models using four measured abundances and 14 upper limits. We assess that models with fit residuals χ²< 3 agree well with observational constraints, resulting in a set of supernova models spanning a wide range of progenitor masses, mixing parameters and explosion energies. The model parameters are however correlated, and the stronger constraints from our improved spectrum modeling do affect the parameter space of models which reproduce observations. We illustrate the co-dependence of parameters in Fig. 10, compared to results from 1D LTE modeling (adopting the same error bars). Our recommended abundances match a more restricted set of models than the 1D LTE abundances, with better matches to low-mass models. The improved modeling is especially restrictive for models with M ≈ 40 M_⊙, where abundances derived from 1D LTE spectrum modeling match supernova models with explosion energies in the very large range 0.6 to 2.7 B. We find that low-mass models (M< 15 M_⊙) matching observations have experienced essentially no mixing (log f_mix ≤ −2.2), and explosion energies are in the range 0.4 to 0.8 B, strongly correlated with mass. More massive models on the other hand explode with energies in the range 1.3 to 2.2 B, while still characterized by little mixing (log f_mix< −1.6). We illustrate the relation, and how it compares to 1D LTE abundance results by overplotting fitted power laws with a break close to M = 14 M_⊙.

We do not include the abundance of Li in our fit to supernova yields, but note that the model with progenitor mass 60 M_⊙ shown in Fig. 9 produces a large overabundance of Li compared to other metals, such that the ejecta mixed with primodial gas is Li-enriched by 0.4 dex. Other models of intermediate mass (40−60 M_⊙) and explosion energy (1.4−2.2 B) produce even more lithium, increasing the abundance by up to 1.6 dex. As Li is created via alpha capture on ³He, produced either via muon neutrino spallation of ⁴He or proton captures on neutrons formed in the charged current reaction (see Heger & Woosley 2010), the yield depends on the adopted neutrino fluxes.

5. Conclusions

We have analyzed the most iron-deficient star known, SMSS0313-6708, using 3D NLTE spectrum synthesis for six elements. We find generally higher abundances than the previous investigation by Bessell et al. (2015), with a significantly higher upper limit for Fe by 1.0 dex, now [Fe/H] < −6.53 (3σ). The increase is mainly due to a more conservative statistical approach (+ 0.5 dex) and improved NLTE data (+ 0.2 dex) rather than the use of 3D NLTE synthesis over ⟨ 3D ⟩ NLTE (+ 0.1 dex), but, nonetheless, all these effects happen to increase the inferred abundance. Our 3D NLTE synthesis also leads to a higher upper limit for Al by 0.5 dex, now [Al/H] < −5.66 (1σ), as well as higher abundances for Mg and Ca by 0.3 dex, now [Mg/H] = −3.83 ± 0.10 and [Ca/H] = −6.94 ± 0.12.

We would like to stress that LTE analyses of extremely metal-poor giant stars, whether in 1D or 3D, are likely to be significantly biased for the elements Mg, Al, Ca and Fe. These abundances may be underestimated by as much as 1 dex. In particular, our 3D LTE analyses of atomic lines are often worse than the 1D LTE analysis in comparison to the 3D NLTE result, due to cancelation effects in the 1D analysis. Such cancelation has previously been predicted (e.g., Asplund et al. 1999) and later found (e.g., Asplund et al. 2003) in analyses of metal-poor dwarfs, and is similar to the so-called NLTE-masking demonstrated in a solar 1D analysis by Rutten & Kostik (1982). Our results from 1D NLTE analyses are similar to those from 3D NLTE, but we find significant differences for the saturated lines of Mg and Ca, which in 1D modeling are sensitive to the v_mic parameter, and for Fe, where the resonance lines are very sensitive to both NLTE and 3D effects. Additionally, 3D−1D corrections inferred in LTE and NLTE−LTE corrections inferred in 1D cannot simply be combined, as the NLTE−LTE difference in line strength varies across the surface of the 3D model. Unfortunately, NLTE models of neutral species are often limited by the atomic data, in particular transition rates due to collisions with neutral hydrogen atoms. These data are now available for a handful of elements (see Barklem 2016b), the bulk of which we utilize in this work.

Combined with the 3D LTE analysis of molecular lines previously reported by Bessell et al. (2015), our 3D NLTE analysis of atomic lines results in abundances being accurately known for all five elements detected in SMSS0313-6708(Li, C, O, Mg and Ca) along with accurate upper limits for four more (N, Na, Al and Fe). Comparing these abundances to predicted supernova yields (Heger & Woosley 2010), we find good agreement with models of progenitor mass 10 M_⊙ and explosion energy <1 B, or 20−60 M_⊙ and explosion energies in the range 1−2 B, where ejecta from the inner core are unable to escape but fall back into the nascent black hole remnant. Explosion energies above 2.5 B can be strongly excluded, as such models eject iron-peak elements producing light-to-heavy element ratios in disagreement with observations.

It is important to extend this work to other ultra metal-poor stars. It is likely that the abundances of other giant stars similar to SMSS0313-6708are significantly underestimated. The magnitude of the 3D NLTE abundance corrections for iron decrease with increasing iron abundance, such that the abundance scale compresses toward higher values, with SMSS0313-6708remaining the most iron-deficient star known by at least an order of magnitude (Frebel & Norris 2015). Correctly inferring the abundance patterns of ultra metal-poor stars, using 3D NLTE modeling, is crucial to constrain the true nature of Population III supernovae, their progenitor stars, and their influence on the early evolution of galaxies.

Acknowledgments

We gratefully thank our referee Matthias Steffen for helpful comments which improved the manuscript, Kjell Eriksson for computing the MARCS model, Alexander Heger for discussions on the supernova pollution scenario and details on the supernova yield grid, and Conrad Chan for help with the starfitcode. T.N. acknowledges support from the Swedish National Space Board (Rymdstyrelsen), and hosting at the ANU for part of this work. A.M.A. and M.A. are supported by the Australian Research Council (ARC) grant FL110100012. K.L. acknowledges funds from the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research as well as funds from the Swedish Research Council (Grant nr. 2015-00415_3) and Marie Sklodowska Curie Actions (Cofund Project INCA 600398). P.S.B. acknowledges support from the Royal Swedish Academy of Sciences, the Wenner-Gren Foundation, Goran Gustafssons Stiftelse and the Swedish Research Council. For much of this work PSB was a Royal Swedish Academy of Sciences Research Fellow supported by a grant from the Knut and Alice Wallenberg Foundation. P.S.B. is presently partially supported by the project grant “The New Milky Way” from the Knut and Alice Wallenberg Foundation. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N under project SNIC2015-1-309 and at UPPMAX under project p2013234.

References

All Tables

All Figures

	Fig. 1 Temperature stratification in the hydrodynamical model shown in a gray-scale histogram as a function of the continuum optical depth. Also shown are its spatial and temporal average (pink line) and a hydrostatic (1D) MARCS model with similar stellar parameters (blue line).
In the text

	Fig. 2 Term diagrams of the atoms analyzed in this work, illustrating all bound-bound transitions considered in the solution of the statistical equilibrium. Thicker lines are due to overlapping transitions. Lines or multiplets used in the abundance analysis are shown in red and labeled. The excited states of Mg ii, Al ii and Fe ii are omitted from the figure.
In the text

Fig. 3 Spatially resolved NLTE equivalent width correction at disk-center intensity (μ = 1), shown for representative lines of Li i (6707 Å, top left), Al i (3961 Å, top right), Ca ii (3933 Å, bottom left), and Fe i (3719 Å, bottom right) in a high-resolution snapshot of the hydrodynamical model (2402 grid points for Li, Al and Ca; 1202 for Fe). Although the behavior is shown for vertical rays, we emphasize that the statistical equilibrium was solved using the 3D radiation field. The NLTE line strength correction, Δlog Wλ,μ = 1 = log Wλ,μ = 1(NLTE)−log Wλ,μ = 1(LTE), corresponds roughly to the abundance correction but with opposite sign. Red and blue colors indicate weaker and stronger lines under NLTE than under LTE, meaning positive and negative abundance corrections. For this illustration, we adopt the abundances A(Li) = 0.70, A(Al) = 1.15, A(Ca) = −0.62, and A(Fe) = 0.60. Na, not shown, behaves very similarly to Li, and the corresponding illustration for Mg is shown in the bottom left of Fig. 4.

In the text

Fig. 4 Spatially resolved line equivalent widths of the Mg i 5183 Å  line at disk-center intensity (μ = 1) in a high-resolution (2402) snapshot of the hydrodynamical model, assuming A(Mg) = 3.80. Top left: NLTE (right: LTE) equivalent widths. Bottom left: NLTE equivalent width correction. Red and blue colors carry the same meaning as in Fig. 3. Bottom right: granulation pattern imprinted on the relative continuum intensity.

In the text

Fig. 5 Variation of predicted disk-integrated line strengths, as inferred from a snapshot from the hydrodynamical simulation downsampled from the original 2402 grid points as well as the horizontally averaged ⟨ 3D ⟩ and hydrostatic 1D case. A vertical dashed line separates the 3D results from the 1D hydrostatic and ⟨ 3D ⟩ results. Calculations in NLTE (black diamonds) and LTE (blue squares) are illustrated for representative lines using A(Li) = 0.70, A(Na) = 0.89, A(Mg) = 3.80, A(Al) = 1.15, A(Ca) = −0.62, and A(Fe) = 0.60. When additional snapshots have been computed (at resolution 602), their individual results are indicated by horizontal lines. A red greek cross indicates the corresponding NLTE result computed in the so-called 1.5D approximation (at resolution 602). Each panel spans over 0.35 dex, except for Al and Fe due to their very strong NLTE effects. For clarity, the line strength determined from the highest-resolution model is indicated as a horizontal dashed line across each panel, and representative errors of 0.03, 0.06 and 0.09 dex are indicated by gray shaded regions.

In the text

	Fig. 6 Upper limit determination for iron, comparing stacked spectra centered on the six strongest unblended lines, Fe i 3440.6, 3581.2, 3719.9, 3737.1, 3820.4 and 3859.9 Å. The 3D NLTE line profiles are shown for the 1σ, 2σ and 3σ upper limits, corresponding to line strengths of 0.7, 1.5 and 2.2 mÅ. The small bump at + 0.4 Å in the synthetic spectrum is due to the nearby line Fe i 3441.0 Å.
In the text

Fig. 7 Abundance results for five lines of magnesium using different types of modeling, shown as a function of the reduced equivalent width. Line central wavelengths are indicated in units of Å. For the 3D NLTE and 3D LTE cases, we overplot linear fits to the trends. In 3D NLTE, we find a positive slope of less than 2σ significance, while 3D LTE as well as the ⟨ 3D ⟩ results indicate >3σ significant negative slopes.

In the text

Fig. 8 Best fitting core-collapse supernova model yields (Heger & Woosley Heger & Woosley 2010) mixed with primordial BBN-composition gas (see text), compared to the inferred abundances and upper limits listed in Table 4. Labeled elements are based on 3D NLTE analysis of atomic lines (blue characters, this work) or 3D LTE analysis of molecular lines (orange characters, from Bessell et al. 2015). Other upper limits are based on 1D LTE analyses (from Bessell et al. 2015). The progenitor stellar mass, explosion energy, and amount of mixing of the best-fit model are indicated in the legend.

In the text

Fig. 9 Best fitting supernova model yields assuming a selection of different progenitor masses, where thicker lines and symbols represent higher mass. The explosion energy (varying color) and mixing fraction (varying plot symbol) have been selected as the best fit to the observed abundances at each selected mass. Residuals are shown beneath for all available elements, with the 1σ, 2σ and 3σ measurement uncertainties indicated in dark to light shades of gray.

In the text

Fig. 10 Comparison of progenitor mass and explosion energy among the well fitted (χ2< 3) supernova models. Colors (explosion energy) and symbols (mixing) carry the same meaning as in Fig. 9. Results are shown on the left for our recommended abundances (from Table 4) and on the right using 1D LTE abundances. To guide the eye, power-laws fitted to the low-mass and high-mass models in the left panel are shown also to the right.

In the text