Solar oxygen abundance

Motivated by the controversy over the surface metallicity of the Sun, we present a re-analysis of the solar photospheric oxygen (O) abundance. New atomic models of O and Ni are used to perform Non-Local Thermodynamic Equilibrium (NLTE) calculations with 1D hydrostatic (MARCS) and 3D hydrodynamical (Stagger and Bifrost) models. The Bifrost 3D MHD simulations are used to quantify the influence of the chromosphere. We compare the 3D NLTE line profiles with new high-resolution, R = 700 000, spatially-resolved spectra of the Sun obtained using the IAG FTS instrument. We find that the O I lines at 777 nm yield the abundance of log A(O) = 8.74 +/- 0.03 dex, which depends on the choice of the H-impact collisional data and oscillator strengths. The forbidden [O I] line at 630 nm is less model-dependent, as it forms nearly in LTE and is only weakly sensitive to convection. However, the oscillator strength for this transition is more uncertain than for the 777 nm lines. Modelled in 3D NLTE with the Ni I blend, the 630 nm line yields an abundance of log A(O) = 8.77 +/- 0.05 dex. We compare our results with previous estimates in the literature and draw a conclusion on the most likely value of the solar photospheric O abundance, which we estimate at log A(O) = 8.75 +/- 0.03 dex.


INTRODUCTION
Oxygen is (behind H and He) the most abundant chemical element in the Universe and it is of major relevance in modern astrophysics, across difference fields, including precision stellar physics, extragalactic astronomy, planet formation, and galaxy evolution. Most importantly, oxygen determines much of the opacity in the solar interior (e.g. Bahcall et al. 2005;Serenelli et al. 2009;Pinsonneault & Delahaye 2009), therefore its abundance is critical to the calculation of Standard Solar Models (SSM), which describe the evolution of the Sun from the pre-main-sequence to the present age of 4.5 Gyr. Oxygen is also the key element in gas-phase spectroscopic E-mail: bergemann@mpia.de diagnostics on extragalactic scales, in particular to infer the metallicities from the H II regions (e.g. Kewley & Ellison 2008;Moustakas et al. 2010) -a broadly-used technique that has been applied to many star-forming galaxies to establish their metal-content and to quantify the mass-metallicity relationship across the entire mass range of galaxies in the local universe. Oxygen is also one of the two -in addition to Mg -most common tracers of nucleosynthesis in massive stars, which explode as core-collapse supernovae. Therefore, O abundance is traditionally used in combination with Fe to map the chemical enrichment and star formation history of stellar populations in the Milky Way and its satellite galaxies (e.g. Tolstoy et al. 2009; Barbuy et al. 2018).
This study focuses on the first problem -the chemical abundance of oxygen in the Sun. Over the past decades, several groups approached this problem from various angles. Some groups attempted to measure the O abundances from the solar photospheric spectrum (e.g. Grevesse & Sauval 1998;Allende Prieto et al. 2001;Asplund et al. 2004;Caffau et al. 2008;Pereira et al. 2009b;Sitnova et al. 2013;Caffau et al. 2015;Socas-Navarro 2015;Steffen et al. 2015;Amarsi et al. 2018). Other groups derived the solar O abundances from the high-ionised lines in the solar wind (e.g. Bochsler 2007;Laming et al. 2017). The latter can be, however, measured less accurately compared to the photospheric analysis methods. In the seminal paper, Grevesse & Sauval (1998) presented the O abundance of 8.83±0.06 dex, based on 1D LTE methods. The most recent estimate, derived by means of a detailed NLTE analysis with 3D radiation-hydrodynamics (RHD) simulations of solar convection, is log A(O) = 8.69 ± 0.03 dex (Amarsi et al. 2018, based on the 777 nm lines), which compares well within the uncertainties with the estimate by Caffau et al. (2015, based on the forbidden 630 nm line) (log A(O) = 8.73±0.02±0.05 dex). Neither of these estimates, however, are satisfactory for the SSM calculations (Villante & Serenelli 2020), as the internal structure of the present-day solar model -the sound speed profile, the depth of the convective envelope -does not compare well with independent constraints on its structure obtained by means of helioseismology. Various scenarios have been put forward to explain the mismatch, ranging from the under-estimated opacity in the interior (Bailey et al. 2015, but see Nagayama et al. 2019) to energy transport by dark matter particles (Vincent et al. 2015). The problem has, so far, not been convincingly explained by any of these scenarios.
In this work, we present a new spectroscopic analysis of the solar photospheric O abundance. We are motivated by the availability of new atomic data and recent 3D magnetohydrodynamic simulations (MHD) of near-surface convection, the chromosphere, and the photosphere (Carlsson et al. 2016). Our model atom of O relies on new data describing collisional excitation and charge transfer with H atoms, photoionisation, and collisions between O i and free electrons. We also analyse the uncertainty in the oscillator strengths of the transitions observed. In addition, we develop a new model atom of Ni in order to study the influence of the Ni blend in one of the diagnostic O i features. Finally, we investigate the O i line formation using different types of the solar model atmospheres, including the Stagger RHD model, but also two 3D MHD models computed self-consistently, with and without chromosphere, using the Bifrost code. This is important for any study that requires precision abundance diagnostics, such as, e.g., a detailed analysis of the Sun relative to solar twins, comparative studies of exoplanet hosts, or spectroscopic age indicators (e.g. Bedell et al. 2014;Buchhave & Latham 2015;Bedell et al. 2018;Lorenzo-Oliveira et al. 2018;Nissen et al. 2020).
The paper is organised as follows. Section 2 provides the details of observational datasets used in the abundance calculations. In Sect. 3 we describe the main properties of the atomic models of O (Sect. 3.1) and Ni (Sect. 3.4), in particular, the calculations of new photo-ionisation cross-sections and electron collision rates (Sect. 3.2), collisions with hydrogen atoms (Sect. 3.3). In Sect. 3.5 we review the 1D and 3D (M)RHD model atmospheres, the statistical equilibrium and radiation transfer codes (Sect. 3.6), and the details of abundance analysis (Sect. 3.7). We then describe the results of our calculations in Sect. 4, addressing the following aspects: line formation of O (Sect. 4.1) and Ni (Sect. 4.3), center-tolimb (CLV) variation (Sect. 4.4), error analysis (Sect. 4.7), the solar O abundance (Sect. 4.6). We perform a comparative analysis of our findings in the context of earlier studies in the literature (Sect. 5). Finally Sect. 6 summarizes the conclusions and provides an outlook for future work.

OBSERVATIONS
In this work, we make use of different spectroscopic data available for the Sun.
Our primary dataset are the new spatially-resolved solar data obtained with the Fourier Transform Spectrograph mounted on the Vacuum Vertical Telescope at the Institut für Astrophysik Göttingen (hereafter, referred to as IAG data). Light from the 50cm-siderostat is picked up by a fiber with an on-sky diameter of 25" and sent to the FTS (for more information, see Reiners et al. (2016) and Schäfer et al. (2020)). The FTS was operated in double-sided mode with a spectral resolution of 0.024 cm −1 , i.e., R = ∆λ λ ≈ 700 000 at λ = 6 000 Å). We collected 123 individual observations of the quiet Sun at 12 different µ-angles during April and August 2020. Each observation took approximately 10 minutes. We used the HITRAN database to identify telluric lines from H2O and O2 and mask them in the spectra. Relative motion between the VVT and the respective solar position was determined using the ephemeris code (Doerr 2015) based on the NASA's Navigation and Ancillary Information Facility SPICE toolkit (Acton 1996). We used the differential rotation law of the Sun (Snodgrass & Ulrich 1990) to correct the individual spectra for solar rotation, and we subtracted barycentric and rotational motion from every observed spectrum. Before addition, we normalized individual spectra and corrected remaining radial velocity offsets between different exposures taken at the same µ-angles but at different positions on the Sun. We verified that these offsets were within the range of pointing errors. Eventually, we added spectra observed at the same µ-angles. Our co-added spectra cover the spectral range between 4700Åand 7900Åand have a signal-tonoise (S/N) of 450 per resolution element at 6000Å. In our analysis we employ the µ angles 1 : 0.2, 0.4, 0.6, 0.8, 1.0.
We also employ high-resolution spatially resolved spectra (Pereira et al. 2009a,b) obtained with the TRIPPEL (TRI-Port Polarimetric Echelle-Littrow spectrograph) spectrograph mounted on the Swedish 1-m Solar Telescope (hereafter, referred to as SST data) (Scharmer et al. 2003). These data have a resolving power of R = λ/δλ ≈ 200 000 and cover the wavelength regions around the O i triplet (7771 to 7775 Å) and [O I] (6300 Å) lines.
In addition, we analyse the spatially resolved solar spectra obtained with the Solar Optical Telescope on the Hinode satellite . These data are not affected by telluric absorption, which is especially critical in the region around 630 nm. The HINODE data have a resolving power of R = 200 000 and a S/N of 20 000. The pointings include the following µ angles: 0.38, 0.70, 0.86, 0.95, and 0.99. For a more detailed description of the observations and the data reduction, we refer the reader to Caffau et al. (2015).
We complement the afore-mentioned data with the very high-quality, R ≈ 350 000, disc-center spectra obtained using the KPNO FTS instrument. These data were released by Brault (1972) and a re-reduced version was published in Neckel & Labs (1984).
The last resource of observational data in this work is the data collected by Takeda & UeNo (2019) using the 60 cm Domeless Solar Telescope (hereafter, referred to as DST data) with the Horizontal Spectrograph at Hida Observatory of Kyoto University. In total, 31 different µ-angles have been observed, of which we were only interested in the disc-center intensity. The spectra have a resolving power of R ∼ 140 000 and S/N of several hundreds.

O model atom
The model atom consists of 122 energy states, with neutral oxygen represented by 120 energy states and the singlyionised oxygen by 2 states, the ground state and the first excited state. The energy levels, which includes levels with principal quantum number n ≤ 30 and orbital angular momentum l ≤ 4, were assembled from the NIST database (Kramida et al. 2020). The original number of O i energy levels is 227, however, we merge the fine structure levels with the energy above 102 700 cm −1 (12.73 eV), so that all states above the level 3s 1 D o are represented by a term with the statistical weight g obtained by summing the individual statistical weights of the levels and the energy computed by weighted averaging of fine structure level energies according to their gvalue (cf. Bergemann et al. 2012). As a result our O i model contains 35 fine structure states and 87 terms. Fig. 1 (top panel) shows the Grotrian diagram of O i, with the key diagnostic transitions in the solar spectrum indicated with bold lines. The majority of lines in the O i system connect energy levels with a very high , E > 10 eV, excitation potential (relative to the ground state).
The diagnostic O lines used in the abundance analysis are the triplet O i lines at 777 nm and the forbidden [O I] line at 630 nm (Table 1). The former lines arise in electric dipole transitions from the 2p 3 3s 5 S o 2 level at 9.15 eV above the ground level. The line at 630 nm is a magnetic dipole transition that arises from the ground level 2p 4 3 P2. This transition is also observed in emission in the airglow spectra in the terrestrial atmosphere (Rees 1989;Slanger et al. 2011). Owing to its very low transition probability (log gf = −9.72, Storey & Zeippen 2000), the [O I] line is weak even in the solar spectrum and its equivalent width (EW) measures only a few mÅ (see Appendix). According to the NIST database 2 , the f-values of the 777 nm and the 630 nm transitions have an uncertainty rating "A" (≤ 3%, or better than 0.01 dex) and "B+" (≤ 7%, or better than 0.03 dex), respectively. However, these uncertainty estimates from NIST appear to be highly overoptimistic. The true uncertainties are likely significantly larger, as indicated by the recently published new transition probabilities Civiš et al. (2018) and our own calculations (Sect. 4.7). In this work, we have chosen to perform all abundance calculations using the average of Hibbert et al. (1991) and Civiš et al. (2018) f-values for the 777 nm lines (log gf = 0.343, 0.197, −0.025 for the 7771, 7774, and 7775 Å lines, respectively), and Storey & Zeippen (2000) data for the 630 nm lines. We furthermore use one half of the difference between Hibbert et al. (1991) and Civiš et al. (2018) (0.026 dex) as a 1σ uncertainty, and adopt a 20% uncertainty (0.08 dex) on the f-value of the magnetic dipole line. We return to this issue in Sect. 4.7. The parameters of the diagnostic lines of O i are provided in Table 1.
The parameters of all other radiative bound-bound transitions in the O model atom, including wavelength, f-values, damping constants, were taken from the Kurucz 3 database. There are 13 483 radiative transitions in the wavelength range from 666 Å to 100 000 Å. However, after merging the energy states and making a cut on the oscillator strength (log gf = −10), the number of radiative transitions is compressed to 1364. The chosen cut is needed to capture the resolved photoionisation and recombination resonances, including the Rydberg Enhanced Recombination (Nemer et al. 2019). These processes lead to a sequence of transitions that involve the diagnostic O I lines, hence they are necessary for the accuracy of the SE calculations. The lines long-ward of 20 000 Å are important, because they represent a channel through which the majority of very high-excitation states are connected and therefore, ensure their convergence to statistical equilibrium (SE).
Most of the lines in the atomic model are represented by Voigt profiles with 9 frequency points, which is sufficient for the SE calculations. The diagnostic O i lines are represented by 251 frequency points. Van der Waals damping, that is, broadening caused by elastic collisions with H atoms, is included using the data from Barklem et al. (2000). It should be noted, though that only 33 of our lines are present in the Barklem et al. (2000) database, therefore for the remainder of transitions we resort to the standard Unsöld formalism (Unsöld 1955). It is known that the Unsöld theory underestimates line broadening, as it only applies to collisional broadening at large atomic separations. However, this is not of a concern for O i, as our tests with Unsöld damping constants scaled by a factor or 2 and 0.5 yield virtually identical results for the SE and the profiles of diagnostic lines.

Photoionisation and electron collision rates
The standard source of photoionisation cross-sections is the Opacity Project (OP) database TOPbase 4 (Cunto & Mendoza 1992). However, the data were computed at a very coarse energy resolution assuming the LS coupling and they are available only for the energy states with the principal quantum number of n ≤ 10. These shortcomings limit the applicability of the OP data to NLTE stellar atmosphere problems.
Also the data that describe the rates of transitions caused by inelastic collisions with free electrons are very sparse. The most recent source of this information is Barklem (2007), who tabulated the rate coefficients for transitions between the 7 low-energy states in O i. They employed the standard R-matrix method assuming LS coupling, which is known to deliver relatively accurate results consistent with direct experimental values. These sparse datasets are, however, not sufficient to describe the entire micro-physics of interactions between O i atoms, radiation field, and electrons. This is especially important, because the key diagnostic O i lines (the 777 nm triplet) connect very high-excitation energy levels ( Fig. 1). Therefore, we opted to compute new photoionisation cross-sections and cross-sections for the transitions caused by e− collisions.
We give a preference to the new photoionisation crosssections over the TOPbase values. First, our target and close coupling expansions are larger than those of the OP. Second, the present cross-sections are in intermediate coupling and account for relativistic corrections, as opposed to the LS-term OP cross sections. Finally, we retained level-to-level partial cross-sections, which allow for a reliable determination of recombination rates. Also, our cross sections have the energy resolution of 2.2 × 10 −4 Ryd, (about 100× better resolution than OP cross sections), which adequately resolves the cross-section resonances. Fig. 2 compares our new values for two selected energy levels with the TOPbase data. Our cross-sections include more resonance structures, owing to the much larger close-coupling expansion. In the case of the 3s 5 S o state, the TOPbase cross-section exhibits an abnormal decreasing trend at higher energies and a sort of bulge at lower energies. Our data show a correct high energy behaviour and reveal that the bulge is, in fact, a series of resonances, which only show up as proper close-coupling channels are included in the calculations. Regarding the cross-section for the 3p 5 P1 level, the results in TOPbase and our data show very similar qualitative behavior, but our cross-section accounts for more series of resonances. It should be noted that TOPbase only provides data in LS coupling, thus the cross-section show in Fig. 2 is actually the state's cross-section.
We also calculated new electron impact excitation rate coefficients using the AUTOSTRUCTURE code (Badnell 2011), which is based on the Breit Pauli Distorted Wave method. All configurations of the form 2s 2 2p 4 , 2s 2 2p 3 nl, 2s2p 5 , and 2s2p 4 nl, with n ≤ 30 and 0 ≤ l ≤ 3 were included in the calculations. For the rest of the O i system, we complemented the R-matrix collision strengths with the rate coefficient computed using the formulae from van Regemorter (1962).

Collisions with hydrogen atoms
Inelastic collisions with H atoms represent an important ingredient in NLTE calculations, as they lead to excitation, de-excitation, ion-pair formation, and mutual neutralization reactions. Detailed quantum-mechanical calculations of the rate coefficients for O were recently presented by two independent groups, Barklem (2018) and Belyaev et al. (2019). These studies rely on different approaches. Barklem (2018) estimated the long-range electronic structure of the OH molecule, by approximating the molecular wave function by the linear combination of atomic orbitals (LCAO) of the atoms that form the molecule, while the calculations of collisional dynamics were performed in the framework of the multichannel model (Belyaev 1993), we refer to this approach as LCAO multichannel study. The O − + H + scattering channel was not taken into account. Independently, Mitrushchenkov et al. (2019) calculated the electronic structure of the OH molecule using the multi-reference configuration interaction (MRCI) method and estimated the rate coefficients for O + H, O + + H − and O − + H + collisional processes using the multichannel model. Later on, Belyaev et al. (2019), based on the MRCI electronic structure by Mitrushchenkov et al. (2019), calculated the rate coefficients for inelastic processes in the same transitions using the quantum hopping probability current (QPC) method. The MRCI QPC study gives a more accurate description of inelastic collisions than the LCAO mul-  Figure 2. New photoionisation cross-sections for the 3s 5 S o 2 and 3p 5 P 1 levels of O i compared with the Opacity Project data (here we take the cross-section in OP LS-coupling of the 5 P for the 5 P 1 level).
tichannel study, because it uses a more accurate electronic structure and take into account both the long-and shortrange non-adiabatic regions, see Belyaev et al. (2019) for the detailed comparison.  Barklem (2018) tabulate the rate coefficients for the transitions between the lowest 18 terms in O i (up to 3s 1 D o at 102584 cm −1 ), whereas Belyaev et al. (2019) provide the data for the lowest 11 terms (up to 4s 3 S o at 96147 cm −1 ). In this work, we complement the latter data with new calculations for 5 additional states of O i with excitation energies from 12.0786 eV to 12.5402 eV, as well as two ionic terms O + + H − and O − + H + . All energy states included in the dynamical calculations are provided in Table 2. For these 5 additional states the transition probabilities are computed using the asymptotic approach (Belyaev 2013), including the non-adiabatic nuclear dynamical treatment accomplished by the multi-channel approach. The calculations are performed within the 4 Σ − molecular symmetry. Higher-lying covalent states are not included in the calculations, because they either create ionic-covalent avoided crossings at inter-nuclear distances greater than ∼ 100 atomic units, and hence the corresponding rates are negligible, or they do not create avoided crossings at all. It should be noted that non-adiabatic tran-sitions between the higher-lying states are possible at shortrange distances, although short-range transitions do not lead to high rate coefficients. These rates can be roughly estimated by means of the Kaulakys analytical approach (Kaulakys 1985(Kaulakys , 1991. These datasets were recently computed by P. Barklem 5 , however, only for the lowest 18 energy states of O i that overlap with the states included in the LCAO multichannel calculations. None of the available datasets accounts for the fine structure, therefore, for the lack of a more accurate formulation, we assign the same rate coefficient to each fine structure level within a given term. We note that an alternative approach was proposed in the literature (e.g. Osorio et al. 2015), however, also this approach involves an assumption of the relative probability of transitions to different final spin states. Our extensive tests suggest that reducing or increasing the quantum-mechanical rate coefficients by a factor of 3 (which corresponds to the average multiplicity of O i) does not have any significant influence on the SE of oxygen, in line with the findings of Bergemann et al. (2019), thus any further arbitrary manipulation of the datasets is not justified.  Barklem (2018), our new estimates, and the values estimated using the model of Kaulakys (1985Kaulakys ( , 1991. Interestingly, there is no systematic difference between the two approaches. For some transitions, Barklem (2018) provide significantly, by up to 4 orders of magnitude, lower rate coefficients. The difference can be explained by the fact that the MRCI QPC study takes non-adiabatic regions at small inter-nuclear distances into account. Due to the transitions at a short range, the probability currents are redistributed differently, as compared to the simplified picture given by including long range non-adiabatic regions (as done in the LCAO multichannel study) only. This also explain the difference in the rates for the important transitions, 3s 5 S o -3s 3 S o and 3s 5 S o -3p 5 P. The rates computed using the Kaulakys (1985Kaulakys ( , 1991 recipe are typically much larger compared to the data computed using the LCAO multichannel or the MRCI QPC studies. The differences are most significant for transitions with ∆E 2.5 eV. In this regime, the Kaulakys model predicts the rate coefficients that are up to 10 orders of magnitude higher compared to the detailed quantum-mechanical calculations. This is not surprising, as the model was developed to describe transitions between high-energy Rydberg states. Fig. 3 (bottom panel) illustrates the rate coefficients of ion-pair formation processes computed using the MRCI QPC and the LCAO multichannel methods. There is no significant difference for the majority of the ion-pair formation processes, except the processes with ∆E < −2.5 eV. But the general behaviour of all ion-pair formation processes is very similar.
It is worth pointing out that the scattering channels j = 9 − 13 belong to the so-called 'optimal window', in accordance with the general behavior found by applying the simplified model (Belyaev & Yakovleva 2017). The model predicts the largest rates for mutual neutralization processes with final-state binding energies in the vicinity of −2 eV, that corresponds to the excitation energy of 11.65 eV for the OH molecule. This optimal window is well understood  (Table 2). Three sets of data, Barklem (2018), Belyaev et al. (2019), and Kaulakys (1991), are shown. through the simplified model and it describes the mechanisms of inelastic atomic collision processes due to long-range nonadiabatic regions created by the ionic-covalent interaction.
Since the LCAO multichannel, MRCI QPC, and Kaulakys data represent physically different approaches to the determination of H impact transition rates, we proceed as follows. We create three NLTE atomic models of O and use them independently in the statistical equilibrium and radiative transfer calculations. One of the models includes the LCAO multichannel data only (hereafter, the LCAO model atom). The other model relies on the LCAO multichannel data co-added with Kaulakys rate coefficients (hereafter, the LCAO+Kaul model), which is appropriate because the LCAO multichannel data describe long-range interactions only and the Kaulakys data are expected to compensate for the lack of short-range interactions. The third model (hereafter, the QPC model) relies on the new rates computed using the MRCI molecular structure and the hopping probability current method. Since in all atomic models, less than 20 energy states (out of 120 O i states in the model atom) are represented with detailed quantum-mechanical data, we connect all other energy states by collision-induced transitions with the rate coefficients computed using the scaled Drawin's formula (Drawin 1968). As we show in Sect. 4.2, the only effect of this approach is that it ensures the collisional coupling with the entire O i system and, therefore, helps to greatly speed up the convergence in the very time-consuming 3D NLTE calculations. Both atomic models also include the rates of charge transfer between the O + and H − ions, computed self-consistently using the methods discussed above.

Ni model atom
Since one of the critical diagnostic O i lines is blended by Ni, we have opted for performing a detailed study of NLTE effects in Ni. To the best of our knowledge, this is the first analysis of the statistical equilibrium of Ni with a comprehensive atomic model with 1D and 3D model atmospheres.
Similar to O, the atomic data for Ni i and Ni ii were adopted from the Kurucz database 6 . The initial datasets comprise 281 energy levels and 9663 radiative transitions for Ni i, and 716 energy levels and 56 193 radiative transitions for Ni ii, respectively. The Ni i levels are taken from the Litzèn et al. (1993) study, which combines experimental measurements with theoretical calculations of the energy level structure. We merge fine structure levels below 6.237 eV and remove transitions with log gf < −10. The uppermost Ni I level in our model has the energy of 7.34 eV, whereas the first ionisation threshold (ground state of Ni II) is located at 7.64 eV. We note that there are only two known energy levels of Ni I that are higher than our threshold for Ni I and are below the 1st ionisation threshold. One of them has the energy of 7.43 eV (59940.517 cm −1 ) and it represents one of the fine structure levels of the 3d8(3F)4s(4F)5p 5 D o term. This level is, therefore, merged with the other lower-energy FS levels of this term. The other level is the 3d8(3F)4s(4F)6s 5 F state with the energy of 7.43 eV (59862.611 cm −1 ). Owing to its small J-value (J = 5) and a very small number of lines connecting it with the rest of Ni I, we merge it with the Ni I state at 7.34 eV. We have verified that our NLTE results do not depend on the treatment of these 2 levels. Our final model of Ni comprises 538 energy states and 2468 spectral lines represented by Voigt profiles with a 9 point frequency quadrature. The first ionisation threshold of Ni i is located at 7.64 eV and that of Ni ii at 18.17 eV. The Grotrian diagram of Ni i is shown in Fig. 1 (bottom panel). Ni has a very complex energy level system with many singlet, triplet, and quintet terms. It should be noted that this diagram only shows energy levels with an LS designation and it is not complete, as many energy levels in Ni i are described in JcK coupling scheme (Litzèn et al. 1993).
Our model atom does not include isotopic splitting and hyperfine structure (for the odd Ni components), as the shifts are very small. Using the data from Johansson et al. (2003), we find, for example, that including the two main isotopes of Ni changes the equivalent width of the Ni i line at 6300.342 Å by less than 0.5%. Therefore, we neglect isotopic shift in the SE calculations of Ni and in the abundance analysis of O.
The photoionisation cross-sections were computed using the hydrogenic approximation. The rates of collisional excitation by electrons and hydrogen atoms were calculated using the standard formulae of van Regemorter (1962) and Drawin (1968). We apply a scaling factor of 0.05 to the rates that describe the collisions of Ni i atoms with H, because the standard formulae are known to severely over-estimate the rate coefficients compared to the detailed quantum-mechanical data (Barklem 2016). However, we investigate in Sect. 4 how this assumption influences the results. The rates of collisional ionisation were computed using the Seaton (1962) formula.

Model atmospheres
We use several types of the solar model atmospheres computed under different physical assumptions. The main difference between the types of models is the treatment of geometry and convective energy transport, as described in the subsequent sections.
For the purposes of comparison with previous calculations in the literature (Asplund et al. 2004;Amarsi et al. 2018), we employ the solar MARCS model also used in Bergemann et al. (2019). This is a 1D LTE line-blanketed model atmosphere computed under the assumption of hydrostatic equi- librium. Convective energy transport is parametrised using the mixing-length-theory (MLT), with the MLT α constant set to 0.5. The micro-turbulence parameter is set to 1 km s −1 .
The 3D model atmospheres employed in this paper were computed using two 3D RHD codes. The STAGGER code (Nordlund & Galsgaard 1995;Nordlund et al. 2009) provides a well established 3D radiation-hydrodynamic simulation of sub-surface stellar convection. The solar simulation used in this work (Collet et al. 2011;Magic et al. 2013) encompasses the entire photosphere as well as the upper part of the convective zone. Bifrost is a 3D (M)RHD code capable of simulating the magnetic solar chromosphere (Gudiksen et al. 2011). We used a Bifrost simulation ranging from ≈ −2.5 Mm below the optical surface up to 8 Mm above. From one snapshot of this full-fledged "chromospheric" simulation, we extracted a purely "photospheric" simulation by cutting off the layers above 1 Mm height. This simulation was given 15 minutes to re-adjust to the new boundary. The horizontal extent of the STAGGER atmosphere is 8x8 Mm, whereas it is 12x12 Mm for the Bifrost atmosphere. Yet both are large enough to host at least 10 granules at the surface, and have the (x,y,z)-resolutions of 240x240x230 (STAGGER) and 512x512x512 (Bifrost), respectively. The horizontal resolution of all models was scaled down to allow for repeated NLTE calculations with modest computational cost. Our detailed tests (Fig. 13, also see Sect. 4.7) allow us to conclude that (x,y) resolutions in excess of (30,30) are not necessary, as already this configuration preserves the main physical properties of line formation and yields line profiles that differ by no more than 2% from the full geometric setup. This uncertainty is sub-dominant to other sources of error in the abundance analysis and can be tolerated. Nonetheless, we account for this error in the determination of the solar O abundance (Sect. 4.6). The vertical resolution remains unaffected, namely 230 grid-points for the STAGGER simulation and 512 or 192 grid-points for the chromospheric and photospheric Bifrost simulations, respectively.
STAGGER and Bifrost codes differ in the input abundances and in the equation of states (EOS). STAGGER is using a modified version of the EOS from Mihalas et al. (1988) and abundances from Asplund et al. (2009). The Bifrost simulation used here was originally made to make comparisons with older simulations from Stein & Nordlund and used an EOS from Gustafsson (1973), abundances and continuum opacities from Gustafsson et al. (1975) computed using the Uppsala background opacity package (Gustafsson 1973).
The temperature structure of the combined down-scaled snapshots of all models is shown in Fig. 4  present the structure of the spatially and temporally averaged ( 3D ) model atmospheres. The averaged models were computed by interpolating the full 3D cubes to an equally spaced optical depth scale and averaging over horizontal slices of equal optical depth. Subsequently all variables were averaged temporally over all spatially averaged snapshots. This was done directly for the temperature and velocity field, but in the case of electron density temporal and horizontal averaging was performed on logarithm of electron number density.

NLTE statistical equilibrium calculations in 1D and 3D
We use two SE codes, MULTI2.3 (Carlsson 1992) and MULTI3D (Leenaarts & Carlsson 2009;Leenaarts et al. 2012), both updated as described in Bergemann et al. (2019) and Gallagher et al. (2020). MULTI2.3 solves the detailed equations of SE in a 1D plane-parallel geometry using the Accelerated Lambda Iteration (ALI) method (Rybicki & Hummer 1991, 1992. Radiation transfer is solved by the method of long characteristics. All calculations in this work were carried out using the local operator acting on the source function. MULTI3D solves the radiative transfer equations in 3D geometry on a Cartesian mesh via the short characteristics method (Kunasz & Auer 1988). It is also capable of solving radiative transfer column-by-column mode, that is ignoring the influence of slanted rays. Extensive tests of the full 3D and the column-by-column solution revealed that the differences between the both approaches are negligibly small, especially for the Sun, as the line EWs change by less than 0.5% (see, e.g. the tests carried out by Bergemann et al. 2019). Our calculation for the O lines also shows that the difference is only 0.6% (Fig. 6). Therefore, we employ the column-by-column radiative transfer solution with MULTI3D in this work. In addition to continuum opacities, both MULTI2.3 and MULTI3D are capable of including background line opacity in the form of a linelist. For O, this not important, because the statistical equilibrium of the element is not sensitive to line blanketing. For Ni i, however, a comprehensive treatment of line blanketing across the entire UV and optical wavelength ranges is critical, because the element, similar to other Fe-group species, is photoionisation-dominated. Since MULTI2.3 is also capable of including sampled opacities, we use the code to compute 3D NLTE corrections for Ni separately using the column-by-column solver and apply them to 3D LTE calculations for Ni, in order to get 3D NLTE line profiles of Ni i lines. This approach is needed, because it is currently not feasible to run MULTI3D with a linelist that is comprehensive and detailed enough, that is, that includes the bound-bound transitions of all relevant absorbers across the entire range of frequencies in the Ni model atom. We have verified (see also Gallagher et al. 2020), that both codes provide identical results when the same initial conditions are used. A detailed analysis of the influence of line blanketing in 3D NLTE calculations for different chemical elements will be presented elsewhere (Semenova et al. in prep).

Abundance calculations for O and Ni
The calculations of O abundance were carried out using the model atmospheres described in Sect. 3.5. In the 3D analysis, we employ 13 consecutive STAGGER snapshots, as well as 10 consecutive snapshots for each of the two Bifrost model atmospheres (with and without chromosphere). All model atmospheres were probed by 17 rays at 5 different angles, one ray being the disc center intensity (cos θ = µ = 1) and the other 16 rays going north, east, south and west at inclinations µ = 0.8, 0.6, 0.4, and 0.2 in order to match the IAG and SST observations.
We use this setup to compute O i line profiles for a dense grid of O abundances, sampling the log A(O) range from 8.50 to 8.90 dex in 3D and 8.50 to 9.20 dex in 1D. Since the diagnostic lines of O are not strong, it is sufficient to assume an equidistant abundance spacing of 0.2 dex. For Ni, we assume the meteoritic abundance of log A(Ni) = 6.23 dex from Lodders (2003). This quantity has an uncertainty of 0.04 dex, and it appears to be superior to the values inferred by spectroscopic methods. Asplund et al. (2009) find the Ni abundance of 6.22 ± 0.04 dex, using 3D LTE modelling. This value is significantly higher than the 3D LTE value recommended by Caffau et al. (2015) (log A(Ni) = 6.11 ± 0.04 dex). The 1D LTE estimate by Wood et al. (2014), based on the analysis of 76 Ni i and Ni ii lines, is 6.27 ± 0.06 dex. We adopt the meteoritic Ni abundance in our work, because of the aforementioned differences in the photospheric results by different groups. We return to this issue in Sect. 5.
The region around the forbidden [O I] line is modelled by co-adding the profiles of the O and Ni lines in 1D LTE, 1D NLTE, or 3D NLTE respectively and fitting the combined line profiles to the observations. Comparing this approximate approach to the exact procedure, in which the absorption coefficients are co-added instead of the intensity or flux profiles, we find that it provides an almost identical solution (differing in EW by less than 0.5%). Our detailed approach to the modeldata comparison and abundance diagnostics is described in Sect. 4.4 and 4.6.

RESULTS
In this section, we describe the results of our calculations using different model atmospheres and model atoms. We start with a brief overview of the formation properties of O i lines in Sect. 4.1, continue with the statistical equilibrium of Ni in Sect. 4.3, discuss our results for the spatially-resolved spectra in Sect. 4.4 and for different atomic models in Sect. 4.5, and comment on the relevance of 3D (M)HD in the abundance calculations in Sect. 4.7.3. We then summarise the methods used in the probabilistic abundance analysis and present the final O abundances in Sect. 4.6. Finally, we compare our results with recent estimates in the literature in Sect. 5 and draw conclusions.

O line formation
The formation of O triplet lines at 7771, 7774, and 7775 Å has been a subject of intense discussions over the past decades. Kiselman (1991) was among the first to point out that for the solar model atmosphere the NLTE effects in the diagnostic O lines, which connect the two high-excitation states, are relatively insensitive to the detailed structure of the model atom and are controlled by scattering in the lines. The formation of O triplet lines in the solar photosphere was extensively discussed in Asplund et al. (2004), Steffen et al. (2015) and Amarsi et al. (2018).
Our calculations confirm that the NLTE effects in the triplet lines are mostly driven by photon losses in the lines (Bergemann & Nordlander 2014) that leads to the deviation of the line source function from the Planck function. The overpopulation of the lower level comes at the expense of the ground state population of O i, via a sequence of transitions that involve charge exchange reactions with O ii, recombinations to the high-excitation O i levels, and spontaneous transitions to the lower levels. This is illustrated in Fig. 7, which shows the departure coefficients bi of the lower 3s 5 S o 2 and upper 3p 5 P energy states, which connect the diagnostic triplet lines, against the continuum optical depth at 5000 Å. Owing to the small Boltzmann factor, the line formation is restricted to a very narrow range of optical depths (−1 < log(τ500) < 0). In this region, the ratio of the departure coefficients of the upper and lower energy states drops below  unity, bj/bi < 1. As a consequence, the NLTE profiles of the permitted O i lines at 777 nm come out stronger compared to the LTE line profiles (Fig. 8 top panel). In contrast, the departure coefficients of the levels involved in the transition (the [O I] line at 630 nm) are very close to unity and the NLTE effects in the line profiles are negligibly small (Fig. 8, bottom panel). This formation of permitted O i lines is very similar in 1D NLTE and in 3D NLTE (see also Asplund et al. 2004;Steffen et al. 2015;Amarsi et al. 2018). The 777 nm lines are stronger in 3D NLTE compared to 1D LTE or 3D LTE, that is, the abundances inferred from 3D NLTE line profiles are significantly lower. The forbidden line is barely sensitive to NLTE effects, even in 3D convective models. This is demonstrated in Fig. 9, which shows the absolute continuum intensities (panel a) and the ratios of the NLTE and LTE equivalent widths computed using the Stagger   Table 3 compares the results in terms of 1D NLTE and 3D NLTE abundance corrections, that is the amount by which the abundance of O has be changed in order to match the EWs of 1D NLTE and 3D NLTE line profiles to those of 1D LTE line profiles. These quantities were computed for the disc-center intensity and fluxes using the O abundance of 8.70 dex 7 . The abundance corrections are not used in the abundance analysis, but are useful for their didactic value. Not surprisingly, the NLTE results and, consequently, the NLTE abundance corrections depend strongly on the input atomic data, foremost on the rates of H impact transitions. Using the quantum-mechanical data (LCAO or QPC), we obtain the 1D NLTE corrections in the range from −0.20 dex (7771 Å line) to −0.16 dex (7775 Å line) for the LCAO model, and −0.17 dex (7771 Å line) to −0.14 dex (7775 Å line) for the QPC model. However, the results change dramatically, if we include the Kaulakys data in the model, by co-adding them with the LCAO datasets. In this case, the 1D NLTE corrections do not exceed −0.11 dex for the 7771 line and −0.08 dex for the 7775 Å line. All other quantities in the atomic models, such as the representation of H-impact transition rates for the high-excitation energy states, electron collisions, or photo-ionisation, do not influence the results at any significant level. Likewise, the NLTE level populations and the pro- 7 We note that the correction is almost identical, if the abundance of 8.90 dex is adopted. files of the diagnostic lines remain nearly identical, regardless of whether Drawin's rates or a blanket constant rate coefficient of 10 −20 are assumed for the majority of the uppermost states in O i (Fig. A1 in Appendix). As emphasised, however, it is important to include these values to ensure rapid convergence, which is important especially in extremely timeconsuming 3D NLTE calculations.

O NLTE abundance corrections
The NLTE corrections computed using the spatially-and temporarily-averaged <3D> Stagger model atmosphere are less extreme. For the LCAO and QPC atomic models (no Kaulakys data), the largest <3D> NLTE corrections do not exceed −0.16 dex and −0.13 dex (the 7771 Å line), respectively. The model atom that includes LCAO and Kaulakys data returns the results that are not too different from LTE: the largest NLTE correction amounts to −0.06 dex (7771 Å line), whereas the 7775 Å line forms nearly in LTE, with the NLTE correction of only −0.02 dex.
Interestingly, the 3D NLTE abundance corrections are very close to the <3D> NLTE results, although the former are slightly larger in absolute value. For the LCAO model atom, the 3D NLTE corrections amount to −0.17 dex for the 7771 Å line, −0.15 dex for the 7774 Å line, and −0.14 dex for the 7775 Å line. The latter line, which is the weakest of all three, is most sensitive to the structure of the model atmospheres and to the properties of the model atom. Its NLTE correction changes from −0.16 to −0.02, depending on the details of modelling. The forbidden oxygen line forms almost in LTE, but it shows a non-negligible sensitivity to convection. The 3D NLTE correction for the 630 nm [O I] line amounts to +0.05 dex relative to 1D LTE. Table 3 also shows our 1D NLTE and 3D NLTE corrections for fluxes computed using the LCAO model atom. These results can be directly compared to the corresponding values from Asplund et al. (2004, their Tab. 2). In 1D NLTE, they obtain −0.24 dex for the 7771 Å line, −0.23 dex for the 7774 Å line, and −0.20 dex for the 7775 Å line. In 3D NLTE, they derive −0.27 dex for the 7771 Å line, −0.24 dex for the 7774 Å line, and −0.20 dex for the 7775 Å line. Our values are in a good agreement with this study, although the amplitude of our NLTE corrections is slightly larger. This can be explained by the complexity of the atomic model. Our model includes 120 states of O i coupled by a variety of processes (Sect. 3.1, whereas the atomic model employed by Asplund et al. (2004) includes only 20 O i states coupled by 43 radiative b-b transitions and electron collisions, and it neglects Himpact collisions. Since the model atom from that study was not published, to check the influence of model atom completeness we performed a sequence of calculations in 1D NLTE (using the MARCS and <3D> model atmosphere), progressively reducing the complexity of the model atom down to 22 levels coupled by 57 radiative b-b transitions. We find that the model that includes 41 oxygen states, with O i closed by the 4f 3 F term (102968 cm −1 ) reproduces exactly the results obtained using the complete model. However, atomic models smaller than this limit lead to more modest NLTE corrections. A 22-level atom underestimates the NLTE corrections by 0.07 dex, resulting into 1D NLTE corrections of the order −0.21 (7771) to −0.16 dex. The residual differences can be explained by Asplund et al. (2004) not including H-impact collisions in their model and using the older version of the Stagger model (M. Asplund priv. comm).
In summary, our findings are similar to those of the pre-

Ni line formation
The statistical equilibrium of Ni was previously investigated by Bruls (1993). Using a model atom with 19 energy levels and 81 radiative transitions, they found that Ni iis subject to over-ionisation and line pumping, the processes that generally favour under-population of energy levels and cause weakening of spectral lines compared to LTE. Over-ionisation and pumping are driven by strong non-local UV radiation field and, similar to other Fe-group species (Bergemann et al. 2012), influence the levels with excitation energies of ∼ 3 to 4 eV. Ding et al. (2002)  Our calculations confirm that Ni i behaves as a classical minority element with numerous important ionisation edges in the UV, which is subject to over-ionisation by strong nonlocal radiation field. As seen in the diagram of the departure coefficients (Fig. 10), the majority of energy levels in Ni i are under-populated. The states with excitation energies of 2 to 4 eV are over-ionised, whereas the higher-lying Ni i states are collisionally coupled with the latter and therefore display similar NLTE effects. The ground state of Ni ii is nearly thermalised across the entire optical depth scale, but higherlying Ni ii states are over-populated, closely resembling the behaviour seen in other Fe-group elements, such as Cr, Ti, Mn, and Co. Table 4 presents the 3D LTE, 1D NLTE, and 3D NLTE corrections for the key Ni line at 6300.34 Å. The line is clearly sensitive to convection, with the 3D LTE estimate being higher by 0.06 dex compared to 1D LTE. The NLTE values were computed using three atomic models of Ni, which differ only in the scaling factor to Drawin's ionisation rates. For the excitation rates, we used the scaling factor of 0.05 (Sect. 3.4). Clearly, assuming no H-impact ionising collisions (S bf = 0) or increasing them by a factor of 100 (S bf = 100) has a significant influence on the NLTE correction, which changes from +0.15 dex to +0.08 dex (MARCS) and from +0.21 dex to +0.14 dex (3D Stagger) respectively. Assuming very large -likely unrealistically too large -Drawin's collisional rates (S bf = 1 000), the difference between 3D NLTE and 3D LTE effectively vanishes. The scaling factors of S bf = 0 and S bf = 1 000 provide the maximum (hereafter, Nimax) and minimum (hereafter, Ni-min) limits on the NLTE effects in Ni lines, respectively. The former approach assumes no H-impact collisional ionisation and the latter results in nearly thermalised ( 5157.98 Å Ni i line, computed using the model atom with collisions scaled by S bf = 1000 and log A(Ni) = 6.23 dex. The latter spectral line was used in the analysis of the solar Ni abundance in Scott et al. (2015) and it appears to be relatively unblended compared to other Ni i features. Both lines are clearly weaker in 3D LTE compared to 1D LTE, and in 3D NLTE the difference is even larger. The weakening of Ni I lines in 3D LTE was also demonstrated by Scott et al. (2015), although no NLTE calculations were performed in that study. To the best of our knowledge, Bruls (1993) and Vieytes & Fontenla (2013) are the only studies of NLTE effects Ni i in the solar spectrum, but they do not provide estimates of line EWs and NLTE abundance corrections in 1D and 3D, so a more quantitative comparison cannot be carried out.
Finally, we note that for Ni i, no quantum-mechanical estimates of photo-ionisation cross-sections, H and e− collision rates are available. So we have to resort to classical recipes, but we caution that these recipes tend to under-estimate the effects of over-ionisation and line pumping, as detailed NLTE studies of other similar species demonstrated (e.g. Mn: Bergemann et al. 2019;Ti: Sitnova et al. 2020). Therefore our present calculations can be viewed as a conservative scenario, yet the actual NLTE effects in the Ni i lines are expected to be larger.

Center-to-limb variation of oxygen abundances
We begin with the analysis of O abundances obtained for different µ angles and for different diagnostic lines of O. The best-fit 3D NLTE and 1D LTE models are compared with the spatially-resolved observations of the Sun in Fig. A2. Figure 12 shows our abundance estimates obtained in 1D LTE, 1D NLTE, and 3D LTE for each of the 5 studied inclinations in the IAG and SST data (Sect. 2). Ideally the slope α of the overdrawn trend line would equal zero, meaning that a given model returns self-consistent results at each angle. Low α values are therefore, to a first order, a good indicator of a more reliable result. Additionally, Fig. A4 (Appendix) shows the observed EWs against the model predictions for 1D LTE, 1D NLTE, and 3D NLTE (LCAO+K). This figure is useful for its pedagogical content, to illustrate the failure of standard 1D or 3D LTE modelling. However, we do not use the actual observed EWs in the abundance analysis (see Sect. 4.6).
Clearly, the 1D LTE assumption does not yield reliable abundance estimates for the O i triplet (Fig. 12, top panels), as there is a strong negative correlation between the individual 1D LTE abundance estimates and the viewing angle, ranging from log A(O) ≈ 8.8 dex at the disc center to log A(O) ≈ 9.3 dex at the limb (µ = 0.2). Also the 3D LTE results are sub-optimal: the slope is slightly smaller compared to 1D LTE, nonetheless the abundances range from log A(O) ≈ 8.85 dex at µ = 1 to ≈ 9.10 dex at µ = 0.2, with a significant systematic bias. 1D NLTE results are similar to 3D LTE in terms of the µlog A(O) slope, however, the disc-center 1D NLTE abundance (here shown for the model atom LCAO+Kaul), log A(O) = 8.66 dex, is in a much better agreement with the 3D NLTE disc-center abundance obtained using the Stagger model, log A(O) = 8.70 dex. However, the 3D NLTE approach clearly outperforms the analyses in that it leads to significantly stronger profiles of the triplet lines at the limb, thereby improving the agreement with the observations (see also Kiselman & Nordlund 1995;Asplund et al. 2004). The results obtained using the IAG and SST data are similar in terms of their CLV abundance slopes.
The forbidden [O I] line (Fig. 12, bottom panels) is less sensitive to the model atmosphere and to the O model atom, as the angle-dependent abundances are similar in 1D LTE and in 1D NLTE. However, all 3D NLTE and 3D LTE values are about 0.10 dex higher compared to 1D results. The differences between various 3D NLTE CLV results are entirely due to the differences in the predicted strength of the Ni blend at 6300.34 Å. The SST results show peculiar undulations, with abundances at µ = 1 and µ = 0.4 being significantly lower compared to other angles. The origin of this systematic is difficult to pin-point. It is possible that the subtraction of fringing and other instrumental artefacts in these data lead to residual systematic biases. In fact, it is clearly seen in the Fig. 2 of Pereira et al. (2009a) that the O I line unfortunately coincides with a minimum in the fringe pattern used to model the continuum of the SST data. The O abundances based on the IAG data at different angles are more internally consistent, supporting the choice of these observations as currently the most reliable dataset.

Influence of O+H collisions
It is useful to look more closely into the abundances computed from different NLTE model atoms of O. In this respect, the most critical reactions are inelastic collisions between O and H atoms, which have a direct influence on the statistical equilibrium of O i.
As shown in Sect. 3.3, our new rate coefficients that describe the transitions in O i caused by collisions with H are different from the previous estimates by Barklem (2018). The latter quantities (LCAO) are on average smaller compared to the data computed using the QPC method. However, if we co-add the data from Barklem (2018) with the data from Kaulakys (1991), in order to account for the lack of shortrang interactions in the LCAO approach (as recommended by Amarsi et al. 2018), the net result is that for the majority of energy levels the combined LCAO+Kaulakys rates are higher compared to the LCAO and to the QPC data.
As a result, the O abundances computed from the 777 triplet lines using the three atomic models show a clear systematic difference (Fig. 12, top panels). The LCAO+Kaulakys model atom leads to a ≈ +0.1 (disc center) to ≈ +0.15 dex (limb) higher solar O abundance compared to the LCAO and QPC models. This is because the latter two models, owing to the less efficient collisional thermalisation they produce, yield larger NLTE effects in the 777 nm O i lines. However, neither of the atomic models can be given a preference, based on the CLV -abundance slope, because the moduli of the slopes are very similar. This suggests that the CLV of O lines alone is not sufficient as a metric to distinguish between the three NLTE atomic models. In other words, each of the model atoms -LCAO, QPC, or LCAO+Kaulakys -allows us to achieve internally consistent results, which are, however, systematically different with respect to each other.
For completeness, we note that the choice of collisional data in the O model atom is of no relevance in modelling the 630 nm oxygen line, as the populations of its energy levels are very close to thermal and the NLTE effects are very minor.

Model-data comparison
We determine the best-fit abundances by employing the χ 2 statistics, comparing the observed data with the model line profiles. Line profiles for a fine grid of abundances are produced and for each observation we find the most probable abundance by minimising the reduced χ 2 : where Npix is the number of wavelength points and N free is 1 as the oxygen abundance is the only free parameter. The observational error σ f obs represents the uncertainty caused by the limited signal-to-noise ratio of the data. For the 777 nm triplet lines, we mask the observations to the region of ±0.5 Å from the line centers and a few surrounding blends (Cr 1 at 7771.76 Å, Ti 1 at 7771.41 Å, Ti 1 at 7772.16 Å, and Sc 1 7772.56 Å), modelled in 1D LTE, are subtracted from the observations. These blends are weak and they are entirely contained in the absorption features associated with the triplet lines. Since the exact position of the line centers varies between different atlases and µ-angles, the models are slightly blue-or red-shifted, so that the synthetic line centers match the observed data.
The 630 nm [O I] line is located in a region, which contains a number of blends, Si 1 (6299.599 Å), Sc 2 (6300.698 Å), and Fe 1 (6301.500 Å). However, models are woefully inadequate to reproduce the observations in this region (Fig.  A5). So to estimate the contribution due to blends, that is to de-blend the O line, we adopted the following approach. We fitted arbitrary Voigt profiles to each of the blending features, allowing for the variation of the amplitude, γ and σ parameters in the Gaussian and Lorentzian profiles. We also allowed for parabolic bisectors in order to deal with the asymmetry of both lines. This fit (Fig. A6) was then used to determine the O abundance from the de-blended O+Ni feature using a mask with the width of ±0.3 Å. The [O I] line was modelled using a grid of oxygen abundances, while the Ni was modelled with a fixed abundance of 6.23 dex. The meteoritic value has an uncertainty of 0.04 dex, but this uncertainty is entirely covered by the difference in the Ni blend contribution caused by adopting different NLTE atomic models of Ni.

Systematic errors
The model (systematic) error σ f model comprises several sources, including the errors caused by the oscillator strengths and damping constants, as well as the errors associated with the limited geometric resolution of the 3D model atmosphere. These errors are not included in the χ 2 calculations (Eq. 1), but are considered separately in the abundance analysis (see Table 5).

Atomic data
The error caused by the uncertainty of oscillator strengths is fully correlated with the abundance error. As described in Sect. 3.1, we use the Hibbert et al. (1991) data for the 777 triplet lines and the Storey & Zeippen (2000) data for the 630 nm forbidden line. However, recently new f values for allowed transitions were presented by Civiš et al. (2018). These data are based on the Quantum Defect Theory (QDT) approach and, for the 777 nm triplet lines they are 12% (−0.05 dex) lower than the f -values provided by Hibbert et al. (1991) and recommended by NIST. The difference between the fvalues from the two sources is four times the NIST estimated uncertainty in the transition.
In order to understand the uncertainties of the f -values we carried out new atomic calculation. We used the code AUTOSTRUCTURE (Badnell 2011) to do multiple calculations of increasing complexity. We started with the same configuration expansion as used by Storey & Zeippen (2000). Further, we added additional configurations with the same orbitals with principal quantum number n≤ 4 and configurations with n=5. In each case we tried different orbital optimization schemes. From the different calculations we find a scatter in the f -value of the 777 nm transitions of about 10% with a tendency for the rate to be closer to lower value of Civiš et al. (2018) than to the Hibbert et al. (1991) value. Therefore, for the 777 lines we adopted the average of both values as our central value, and the associated f -value error (2 σ) is assumed to be the difference between the both quantities.
The f -value for the magnetic dipole line at 630 nm is much more difficult to pin point. Such weak lines are generally very sensitive to cancellation effects among different contributions from configuration interaction representations of the atomic levels. The main challenge is not the convergence of the expansion, but the sensitivity to how the orbitals are optimized. The NIST recommended transition probability for this line is 5.63 × 10 −3 s −1 with the uncertainty flag B+ (better than 7%), which is adopted from Baluja & Zeippen (1988) and Fischer & Saha (1983). The most recent calculations by Storey & Zeippen (2000) yield 6.45×10 −3 s −1 . This is a difference of ∼ 15%, which is more than twice the estimated uncertainty in this rate by NIST. The results of our own calculations show the rate varies widely between ∼ 3 × 10 −3 s −1 and ∼ 1 × 10 −2 s −1 , with the best value being ∼ 8 × 10 −3 s −1 . As stated earlier, we use the Storey & Zeippen (2000) rate in our analysis, however, to account for the extreme sensitivity of the f -values to the details of orbital optimisation, we assume a more realistic uncertainty of 20%.
Finally, we investigated the uncertainty associated with damping, by performing spectrum synthesis calculations with a modified, by ±10%, value of the damping constant. This, however, has a negligible influence on the O triplet and [O I] lines and is therefore neglected.

Geometric resolution
We also take into account a systematic error associated with the limited geometric resolution of 3D model atmospheres used in radiative transfer calculations. This correction is estimated as follows. The two families of photospheric 3D atmospheres employed in this paper have (x,y,z)-resolutions of 240x240x230 (STAGGER) and 512x512x192 (Bifrost-phot). We restrain from re-sampling the vertical resolution, but down-sample the number of vertical columns, ensuring that we have enough columns to maintain the ratio of intergranular lanes to granules covering the entire simulation. We converted a single STAGGER snapshot to horizontal resolutions 5x5, 10x10, 20x20, 30x30 and 120x120 and computed equivalent widths of the diagnostic O i lines for each downsampled snapshot. Similarly, equivalent widths for a photospheric Bifrost snapshot with different horizontal resolutions between 10x10 and 160x160 were computed. These calculations were carried out in LTE, because our previous analysis (Bergemann et al. 2019)  x-direction y-direction Figure 13. Influence of the geometrical resolution (nx,ny) of 3D STAGGER model on the measured equivalent widths. Here nx is kept constant at 5. ny varies from 5 to 240, where 240x240 corresponds to the full resolution of the model. The vertical dashed line shows our default geometric resolution (30,30) adopted in 3D NLTE radiative transfer calculations. resolution of (x,y)= (30, 30) offers a reasonable compromise between the computational expense of 3D NLTE modelling and the physical realism, leading to a small (∼ 1 − 2%) systematic error in the line strength compared to that obtained the full geometric setup. To correct for this small bias, we extract the geometric resolution correction (Table 5) by comparing the EWs of lines computed at our nominal resolution and at the highest-possible resolution. This correction is then separately applied to the abundance determined at each µ angle. The results are shown in Fig. 13. The disc-center and limb results are not very sensitive to this effect, but at µ = 0.8 and µ = 0.6, the EW error caused by using the horizontal resolution of (x,y)= (30, 30) can be as large as 3%. For the disc-center, the abundance correction for the diagnostic O lines amounts to 0.01 dex. We do not use the Bifrost models to derive the absolute solar O abundance (Sect. 3.5), but only to quantify the influence of the chromosphere, that is, the re-adjustment of the photospheric structure under the influence of the chromosphere, on the line formation and on the O abundances. This is done by comparing the results calculated with two types of Bifrost models: (a) the full simulation covering the photosphere and the chromosphere, and (b) the strictly photospheric simulation (Sect. 3.5). Each simulation represents an extended time series, from which we extract 10 snapshots for the detailed radiative transfer. The average difference between the T eff values of the snapshots is 10 K.
Comparing the O abundances derived using the two Bifrost model sequences in 3D NLTE, we find we find that the abundance differences are not large (Table 5). The chromosphere weakens the 777 triplet lines that corresponds to the increase of O abundance by ∼ +0.01 dex at the disc center and +0.02 dex at the limb relative to a pure photospheric simulation. The 6300 [O I] line is relatively unaffected by the chromosphere, and the O abundance inferred from this line is only 0.003 dex lower. We note that even ignoring the snapshot with the largest T eff difference (16 K) does not alter the estimate of chromospheric correction by more than 0.001 dex.
Whereas the chromospheric correction for the O lines is small, the presence of a chromosphere may have a larger influence on the spectral lines that form higher up in the layers close and above the temperature minimum, such as the Ni 6767 Å helioseismology line. Table 6 provides our estimates of the photospheric O abundance computed in 3D NLTE from different observational datasets, where the results were corrected for the chromospheric and spatial resolution effects (Sects. 4.7.2,4.7.3). The individual results obtained in 1D LTE, 1D NLTE, 3D NLTE with different atomic models are provided in Tab. A1 (Appendix). The EWs of the corresponding lines are tabulated in Tab. A2 (Appendix).

Final O abundance
To compute the final estimates of O abundance and its uncertainties from different observational data and different atomic models, we proceed as follows. For both 777 and 630 nm lines, we adopt the central values as derived from the IAG data using the LCAO+K model atom of O, because of the superior resolving power of these observations and more consistent estimates of abundances derived from the spatiallyresolved spectra taken at different pointings across the solar disc 8 . The LCAO+K (with Kaulakys data) atomic model, furthermore, provides a flatter CLV abundance slope, compared to the LCAO or QPC collisional data, but also from the perspective of atomic physics, the Kaulakys data appear to be necessary, as they are expected to compensate for the lack of short-range interactions in the quantum-mechanical calculations. We also adopt the S bf = 100 model atom of Ni, as this NLTE model provides the average of the plausible 3D NLTE Ni blend contributions to the [O I] feature that corresponds to the average between the extremely high, moderate, and negligible collision rates (Sect. 4.3).
To estimate the final combined uncertainties, we resort to a Monte Carlo simulation and make use of the central limit theorem, which states that the normalised sum of random uncorrelated variables closely follows a normal distribution, regardless of the shape of the distribution function of individual variables. This is a commonly used approach in case when the shapes of uncertainties are not known, as in our case. The errors associated with the analysis (see Tab. 5) are used to generate three uniform distributions (large samples of 1 million points each), representing (a) the uncertainty of the oscillator strengths (f -value), (b) the uncertainty caused by the Ni blend (630 nm line) and that caused by the O+H collisional data (777 nm lines), and (c) the uncertainty caused by using different observed spectra. The uncertainty caused by the collisional data is taken to be the difference between the QPC and LCAO+K results. The latter is assumed to be one half of the difference between the SST (min) and the IAG (max) results for the 777 nm line, and that between the SST (min) and DST (max) results for the 630 nm line. These errors are co-added and one standard deviation (1σ) of the resulting distribution (Fig. 14) is adopted as the final uncertainty. We note that the final combined distribution of errors for the 777 nm line is very close to normal, whereas for the 630 nm line there is a small deviation owing to a large dominant uncertainty source. However, this has no influence on the conclusions. As a result, our 3D NLTE value derived from the 777 and 630 nm lines is log A(O) = 8.739 ± 0.027 dex and log A(O) = 8.771 ± 0.053 dex, respectively. Assuming that the results are independent, we can combine them to obtain the final photospheric O abundance value of log A(O) = 8.75 ± 0.03 dex.
We note that if we were to use the old Hibbert et al. (1991) f -values and the SST data, our results based on the LCAO+K model atom of O would be in excellent agreement with the 3D NLTE estimate by Amarsi et al. (2018) who relied on the same input. Likewise, our 1D NLTE result based on the MARCS model is in agreement with that by Sitnova & Mashonkina (2018).

DISCUSSION
Our estimates of O abundances based on the permitted 777 nm triplet lines and on the 630 forbidden [O I] line are in a good agreement within the respective uncertainties. The triplet lines are, however, sensitive to the input physics of the O model atom, whereas the forbidden line is sensitive to the input physics of the Ni model atom, primarily to the collisional data. Also, the f -values carry a significant source of uncertainty, as the direct experimental verification of the transition probabilities for the 777 nm and 630 nm is not possible. These errors are associated with the physical limitations of theory of atomic physics, molecular physics, and collisional dynamics, and, therefore, their true distributions are unknown. Nonetheless, considering all sources of error, it is very plausible that the photospheric O abundance is close to log A(O) = 8.75 dex.
It is interesting to discuss our results in comparison with the other recent estimates of the solar photospheric O abundance. Amarsi et al. (2018) and Asplund et al. (2021) found log A(O) = 8.69 ± 0.03 dex and 8.69 ± 0.04 dex, respectively. The former estimate is based on the 777 nm line, while the latter also includes the 630 nm and selected molecular OH lines. In principle, our results and those by this group are not inconsistent within their combined uncertainties. Our value is in a much better agreement with the values proposed by Caffau et al. (2008), Caffau et al. (2015), and Steffen et al. (2015), who analysed different O lines, including several forbidden [O I] lines, the 777 triplet, and the near-IR O I features. Another recent study of the solar photospheric O abundance is that by Cubas Armas et al. (2020). Their analysis is based on a semi-empirical method of determining the abundance by inverting spatially-resolved observed solar spectra. Their estimate, based on the forbidden [O I] line at 630 nm, is log A(O) = 8.80 ± 0.03 dex, also consistent with our results for this feature.

CONCLUSIONS
We present a re-analysis of the solar photospheric O abundance using new atomic data for O and different 3D radiation hydro-dynamical models of solar atmospheres. The NLTE model atom includes novel rates of collisions between O and H atoms, computed using the quantum hopping probability current method , increasing the number of O i terms with accurate collisional data to 16. The R-matrix method for atomic scattering calculations (Berrington et al. 1978) is employed to compute new photoionisation cross-sections and the rates describing transitions caused by inelastic processes in collisions of O atoms with free electrons. We also present a new comprehensive model atom of Ni, which, for the lack of detailed quantum-mechanical data, is still based on classical formulae for photoionisation and collisional reaction rates. We use the 1D MARCS solar model atmosphere, the 3D Stagger model (Collet et al. 2011;Magic et al. 2013), and two version of 3D MHD models computed self-consistently, with and without the chromosphere, using the Bifrost code (Carlsson et al. 2016).
The model atoms are used to performed 1D LTE, 1D NLTE, and 3D NLTE radiative transfer calculations and abundance analysis of O lines in the solar spectrum. We focus on the least blended lines of O i at 7771, 7774, 7775 Å, and we also consider the forbidden O i line at 6300 Å. We also perform 3D NLTE calculations of the Ni i blend at 6300.341 Å, which contributes about 25% of the entire absorption feature. The observational data are taken from different sources. Our primary source of data are the new spatially-resolved spectra taken with the IAG FTS instrument. The data have a resolving power of R ∼ 700 000, much higher compared to all previously available solar data. Similar to previous studies, we also include the spectra obtained with the ground-based Swedish Solar Telescope (Pereira et al. 2009a), the Domeless Solar Telescope (Takeda & UeNo 2019), and with the space-based Solar Optical Telescope on the HINODE satellite .
We obtain the solar O abundance of log A(O) = 8.739± 0.027 dex and log A(O) = 8.771± 0.053 dex based on the 777 nm triplet and 630 nm forbidden [O I] lines, respectively. Our final combined value is log A(O) = 8.75 ± 0.03 dex, which is close to the estimate proposed in the compilation of the solar abundances by Caffau et al. (2008, 8.76±0.07 dex), but slightly higher than the recent value by Asplund et al. (2021, 8.69 ± 0.04 dex). The differences between our and previous results can be explained by the updated transition probabilities for the 777 lines, new higher-resolution observations collected with the IAG FTS facility, 3D NLTE modelling of the Ni blend, and a correction for a systematic bias caused by the chromospheric back-heating. These effects add up to increase the solar abundance by about +0.05 dex relative to the value proposed by Asplund et al. In addition, we confirm that using the HINODE space-based data, which were employed by Caffau et al. (2015), we recover a slightly higher abundance compared to the value obtained by using the SST data chosen by Asplund et al. (2021).
Our results thus suggest that the question of whether the solar photospheric O abundance is low or high, that is closer to the values from Grevesse & Sauval 1998, is still open. There are remaining controversies associated with modelling of the permitted 777 nm lines and of the forbidden [O I] line. First, the large systematic differences between two sources of Himpact collisional data (Kaulakys 1991 andBelyaev et al. 2019) must be resolved. Second, the differences between the new lower f -values from Civiš et al. (2018) and the older values from Hibbert et al. (1991) must be settled. Until then, we recommend to use the average of the values from both sources. Our results suggest that the indirect influence of the chromosphere on the O lines in the solar spectrum is not large. However, it is desirable to perform new self-consistent 3D R(M)HD simulations of the solar atmosphere, including chromosphere, with state-of-the-art micro-physics and Non-LTE radiative transfer. Figure A1. Profiles of the 7771 Å line computed in NLTE using three different atomic models with the collisional rate coefficients for the upper 100 states uniformly set to 10 −20 , assumed to follow the Drawin's recipe, and scaled by a factor of 10 relative to Drawin.