A&A 456, 675-688 (2006)
DOI: 10.1051/0004-6361:20064986

Line formation in solar granulation

VII. CO lines and the solar C and O isotopic abundances[*]

P. C. Scott1 - M. Asplund1 - N. Grevesse2,3 - A. J. Sauval4

1 - Research School of Astronomy and Astrophysics, Mt. Stromlo Observatory, Cotter Rd., Weston Creek, ACT 2611, Australia
2 - Centre Spatial de Liège, Université de Liège, avenue Pré Aily, 4031 Angleur-Liège, Belgium
3 - Institut d'Astrophysique et de Géophysique, Université de Liège, Allée du 6 août 17, B5C, 4000 Liège, Belgium
4 - Observatoire Royal de Belgique, avenue circulaire 3, 1180 Bruxelles, Belgium

Received 9 February 2006 / Accepted 18 April 2006

CO spectral line formation in the Sun has long been a source of consternation for solar physicists, as have the elemental abundances it seems to imply. We modelled solar CO line formation using a realistic, ab initio, time-dependent 3D radiative-hydrodynamic model atmosphere. Results were compared with space-based observations from the ATMOS space shuttle experiment. We employed weak 12C16O, 13C16O and 12C18O lines from the fundamental ( $\Delta v = 1$) and first overtone ( $\Delta v = 2$) bands to determine the solar carbon abundance, as well as the 12C/13C and 16O/18O isotopic ratios. A weighted solar carbon abundance of $\log\epsilon_{\rm C}=8.39$ $\pm$ 0.05 was found. We note with satisfaction that the derived abundance is identical to our recent 3D determination based on C I, [C I], C2 and CH lines, increasing our confidence in the accuracy of both results. Identical calculations were carried out using 1D models, but only the 3D model was able to produce abundance agreement between different CO lines and the other atomic and molecular diagnostics. Solar 12C/13C and 16O/18O ratios were measured as 86.8 +3.9-3.7 ( $\delta^{13}{\rm C} = 30^{+46}_{-44}$) and 479 +29-28 ( $\delta^{18}{\rm O} = 41^{+67}_{-59}$), respectively. These values may require current theories of solar system formation, such as the CO self-shielding hypothesis, to be revised. Excellent agreement was seen between observed and predicted weak CO line shapes, without invoking micro- or macroturbulence. Agreement breaks down for the strongest CO lines however, which are formed in very high atmospheric layers. Whilst the line asymmetries (bisectors) were reasonably well reproduced, line strengths predicted on the basis of C and O abundances from other diagnostics were weaker than observed. The simplest explanation is that temperatures are overestimated in the highest layers of the 3D simulation. Thus, our analysis supports the presence of a COmosphere above the traditional photospheric temperature minimum, with an average temperature of less than 4000 K. This shortcoming of the 3D model atmosphere is not surprising, given that it was never intended to properly describe such high layers.

Key words: convection - line: profiles - Sun: abundances - Sun: photosphere - Sun: infrared - solar system: formation

1 Introduction

1.1 CO in the solar spectrum

Thousands of vibration-rotation lines from the fundamental $^1\Sigma^+$ electronic band of CO are present in the infrared solar spectrum obtained from space in the Atmospheric Trace Molecule Spectroscopy (ATMOS) space shuttle mission (see Sects. 3 and 4). From 1350 to 2328 cm-1 (i.e. 4.3 to 7.4 $\mu$m; fundamental bands from 1-0 to 20-19) and 3410 to 4360 cm-1 (i.e. 2.3 to 2.9 $\mu$m; first overtone bands from 2-0 to 14-12), the solar spectrum looks more like a pure CO absorption spectrum, displaying not only 12C16O lines but also 13C16O, 12C18O and 12C17O isotopomeric lines. Thanks to their very high excitations (higher than in any laboratory spectra where rotational excitations are concerned), the solar CO lines have been used to derive a new set of highly accurate molecular constants (Farrenq et al. 1991; Sauval et al. 1992). As the most sensitive indicator of solar photospheric temperatures, and formed over an extremely wide domain of optical depths (Fig. 1), analysis of the CO lines offers a unique opportunity to test physical conditions up to very high atmospheric layers, derive isotopic abundance ratios and refine the solar carbon abundance (Grevesse & Sauval 1994,1991; Grevesse et al. 1995; Grevesse & Sauval 1992). Only a limited number of clean CO lines are available in ground-based spectra however, with the $\Delta v = 1$ sequence heavily masked by strong, broad telluric absorption and the $\Delta v = 2$ sequence also strongly polluted by terrestrial lines.

1.2 CO analyses and problems

Since the seminal observations of Noyes & Hall (1972) revealed limb CO core brightness temperatures below the classical photospheric temperature minimum of 4500 K, debate has raged over what CO absorption features in the solar spectrum really tell us. Ayres & Testerman (1981) found that the most natural explanation was an absence of any photospheric temperature minimum, requiring an outwardly decreasing temperature structure from the base of the photosphere right through the chromosphere. However, given that Ca II and Mg II line cores indicate chromospheric temperatures well above the classical minimum (Ayres & Linsky 1976), Ayres (1981) proposed the existence of a cool CO structure in the low chromosphere, coexistent with other areas of hotter gas responsible for the Ca II and Mg II core emission: a "thermally bifurcated'' low chromosphere.

Ayres & Wiedemann (1989) went on to test the LTE assumption of the bifurcation model, finding that non-LTE effects had <2$\%$ effect upon CO line cores, even as near to the solar limb as $\mu = 0.1$. Whilst these findings were disputed (Mauas et al. 1990), the existence of some cool CO in the atmosphere was accepted, though temperature inhomogeneities were suggested to be due to the effects of localised mechanical heating rather than CO radiative cooling. The presence of cool CO in the chromosphere was put beyond doubt four years later with the discovery of off-limb CO emission indicating temperatures of less than 4000 K in some places (Solanki et al. 1994).

Using imaging spectroscopy, Uitenbroek et al. (1994) studied spatial and temporal variations in CO absorption at disk centre, observing temperature inhomogeneities of up to 600 K correlated with hydrodynamic perturbations. This result supported the existence of the temperature bifurcation, though driven more by dynamic processes than CO cooling and magnetic heating as suggested originally. This was picked up on in a subsequent revision of the bifurcation model (Ayres & Rabin 1996) in which the cool component was dubbed the "COmosphere''.

Carlsson & Stein (1995,1997) demonstrated that chromospheric Ca II H2V and K2V grains are formed by acoustic shocks, using a 1D non-LTE hydrodynamic model. Their results strongly supported the notion pushed by Ayres & Rabin (1996) of a bifurcated atmosphere with differences caused by dynamic effects, even also suggesting some influence by CO cooling owing to the inclusion of instantaneous chemical equilibrium (ICE) CO opacities in the model. Kalkofen et al. (1999) and Kalkofen (2001) argued against the Carlsson-Stein model because it did not adequately represent the entire solar energy output, and so could not reproduce the observations of ubiquitous chromospheric emission. However, this was never Carlsson & Stein's intention, as discussed by Ayres (2002) in an extensive rebuttal of Kalkofen et al.'s criticism. Ayres went on to show that despite not being designed with CO diagnostics in mind, the Carlsson-Stein model far outperformed Kalkofen et al.'s preferred model (VAL; Vernazza et al. 1976,1981,1973) when put to work on solar CO lines.

Meanwhile, Uitenbroek (2000b) performed detailed 1D NLTE simulations of CO absorption in the solar disk, reconfirming that the LTE approximation is appropriate in this context. Uitenbroek (2000a) also studied 3D LTE CO line formation using a single snapshot from an early version (Stein & Nordlund 1989) of the current state-of-the-art 3D atmospheric models (Asplund et al. 2000b; Stein & Nordlund 1998). The calculated line profiles were deeper than those observed by the ATMOS mission, something Uitenbroek attributed to the ICE approximation.

Growing concern (Ayres 2002; Uitenbroek 2000a,b) over the appropriateness of the ICE approximation was addressed by Asensio Ramos et al. (2003), who investigated departures from ICE in CO lines using the Carlsson-Stein model. Results showed significant concentrations of cool CO gas extending just 700 km above the Sun's surface, with the ICE assumption valid below this height. 2D simulations with non-equilibrium chemistry (Wedemeyer-Böhm et al. 2005) have confirmed that the ICE approximation is valid in the photosphere and low chromosphere. Despite employing non-equilibrium chemistry in both the hydrodynamic and radiative transfer simulations, modelled CO lines were still considerably deeper than observed profiles. This is probably because the temperature structure of the excerpt from the 2D model was too cool, producing an effective surface temperature more than 300 K lower than the true solar value. Furthermore, 2D simulations tend to give erroneous photospheric line profiles compared with 3D calculations because of the geometric restriction on the convective flow (Asplund et al. 2000a).

1.3 Solar C and O isotopic abundances

The solar carbon abundance has recently been revised (Asplund et al. 2005a) from the $\log\epsilon_{\rm C} = 8.52$ $\pm$ 0.06 given by Grevesse & Sauval (1998). Analyses of [C I], C I, CH and C2 (Allende Prieto et al. 2002; Asplund et al. 2005b) now indicate that the logarithmic solar carbon abundance is in fact 8.39 $\pm$ 0.05. These results provide a comparison for the solar carbon abundance determination we carry out using CO lines. Considering the difficulty of reproducing observed CO lines with theoretical calculations, these provide a very stringent test of the model atmosphere. Asymmetries of the CO lines constitute an even greater challenge for theoretical atmospheres (Blomme et al. 1994).

Early measurements of the solar 12C/13C ratio using CH lines were typically higher than the terrestrial reference ratio of 89.4 $\pm$ 0.2 (Coplen et al. 2002), and exhibited large uncertainties (see Harris et al. 1987 for an inventory of early measures). The most reliable (lowest uncertainty) determinations to date utilised CO lines, producing near-terrestrial ratios of 84 $\pm$ 9 (Hall 1973) and 84 $\pm$ 5 (Harris et al. 1987). The solar 16O/18O ratio was also measured, producing values of >500 and 440 $\pm$ 50 respectively. These are similar to the representative terrestrial ratio of 498.7 $\pm$ 0.1 (Coplen et al. 2002)[*].

Different lunar regolith (surface) analyses of the solar wind have indicated $\delta^{13}{\rm C}\leq -105$ $\pm$ 20 (Hashizume et al. 2004) and $30\leq\delta^{13}{\rm C}\leq-30$ (combined results reviewed in Wiens et al. 2004). $\delta^{18}$O values tabulated in the same review for direct measurements of solar wind particles by two different satellites give 110 +450-250 and 120 +280-190. For comparison, the results of Harris et al. (1987) on the same scale imply $\delta^{13}{\rm C} = 59^{+67}_{-60}$ and $\delta^{18}{\rm O} = 130^{+145}_{-115}$, without the inherent uncertainty in the scale taken into account. The most recent lunar regolith analyses differ in their estimates of solar  $\delta^{18}$O, with one suggesting a negative value (Hashizume & Chaussidon 2005, see also Davis 2005) and the other producing $\delta^{18}{\rm O} = 50$ (Ireland et al. 2006). Fractionation that might occur in the solar wind is still not properly quantified however (Wiens et al. 2004), so such results may say as much about this process as they do about the actual photospheric ratios. Solar chemical evolution is also thought to effect abundance agreement between the current photosphere and lunar inclusions irradiated by the prehistoric solar wind, though the predicted difference over the Sun's lifetime is less than 10 permil in  $\delta^{13}$C or  $\delta^{18}$O (Turcotte & Wimmer-Schweingruber 2002). Recently, Yurimoto & Kuramoto (2004, see also Yin 2004) proposed a model of solar system formation whereby oxygen isotopic fractionation occurs in the protosolar cloud due to UV  irradiation. This theory predicts a solar  $\delta^{18}$O of -50, something that can be directly tested with the 18O/16O measurement to be performed herein.

1.4 This study

In this paper, we investigate CO absorption line formation in the solar spectrum using a three-dimensional radiation-hydrodynamic simulation of the Sun's atmosphere. This analysis is similar to that of Uitenbroek (2000a), but with an improved 3D model as its basis. This enables a better description of CO and a more precise evaluation of the reality of the model in the context of CO, as well as simulation of line formation in the disputed COmosphere. It should be noted that the model is not intended to reproduce the ubiquitous chromospheric emission nor near-limb CO lines, as it does not include magnetic fields nor extend to great enough heights.

Using the same model and data, we also arrive at new estimates of the solar carbon abundance, 12C/13C and 16C/18O ratios using appropriate sets of weak 12C16O, 13C16O and 12C18O lines. Unfortunately, the 12C17O lines clearly identified in the ATMOS spectra are too weak and perturbed to permit any defensible estimate of the solar 16O/17O ratio, in our opinion. The isotopic ratios constitute an improvement over the 1D analyses of Hall (1973) and Harris et al. (1987), thanks to the use of 3D model atmospheres, more recent gf values, more appropriate line lists and better observations. We believe that our ratios are also superior to alternative determinations carried out very recently by Ayres et al. (2006), as will be detailed in Sect. 7. The determination of the solar carbon abundance using CO lines provides an important comparison with the study of Asplund et al. (2005b) based on [C I], C I, CH and C2 lines. It also circumvents past trouble in using CO lines and 1D models to measure the solar carbon abundance (e.g. Grevesse et al. 1995).

The 3D model and spectral line lists analysed are briefly described in Sects. 2 and 3 respectively, whilst Sect. 4 details our manipulation of the ATMOS data. The study of the detailed shapes of strong CO absorption lines and their implications for the temperature structure in very high atmospheric layers is presented in Sect. 5. In Sect. 6 we present abundance and isotopic ratio determinations, based on weak CO lines. Comparisons with earlier work are given in Sect. 7, and our conclusions are summarised in Sect. 8. A detailed description of the apodization procedure applied to modelled spectra is given in Appendix A. Appendix B provides derivations of scaling factors used in the analysis, and Appendix C contains full line lists.

2 Model atmospheres and line formation calculations

The hydrodynamic simulation used was described by Asplund et al. (2000b). Specifically, it covers a physical area 6 $\times$$\times$ 3.8 Mm of which about 1 Mm is above the optical solar surface ( $\tau_{\rm Ross} \approx 1$), at a resolution of 200 $\times$ 200 $\times$ 82. The simulated domain was bounded below by a transmitting boundary, and above by an extended transmitting boundary across which the density gradient was kept hydrostatic. Horizontal boundaries were periodic. The MHD equation-of-state (Hummer & Mihalas 1988; Mihalas et al. 1988) and Uppsala opacities were used, with LTE assumed. Continuous and line opacities were calculated using opacity binning. The model included no tuneable free parameters, and was characterised as solar by the accepted gravity, effective temperature and standard solar composition of Grevesse & Sauval (1998). A numerical viscosity was employed to stabilise the simulation. About 100 snapshots of the convective simulation were stored, representing approximately 50 min of real solar time. Radiative transfer calculations for spectral line formation were carried out over an interpolated 50 $\times$ 50 $\times$ 82 grid, using the Uppsala equation-of-state and opacities with molecular densities determined under the ICE approximation. LTE was again assumed, as validified for the CO lines by Ayres & Wiedemann (1989) and Uitenbroek (2000a). This is the same atmosphere and line formation code recently used to re-derive the solar abundances of all elements between Li and Ca (e.g. Asplund 2000,2004; Asplund et al. 2005b,2004,2000c,2005a; Asplund et al., in preparation; Allende Prieto et al. 2002,2001).

For comparative purposes, abundance calculations were carried out using four different model atmospheres: the 3D hydrodynamic model, the 1D HM model (Holweger & Müller 1974), the 1D MARCS model (Gustafsson et al. 1975; Asplund et al. 1997), and a contraction of the 3D model into the vertical dimension only, which we designate "1DAV''[*]. The horizontal averaging used to produce the 1DAV model was performed over surfaces of common optical depth rather than geometrical height, as the optical depth scale has most relevance to line formation. All three 1D models included microturbulence of 1 km s-1. All three were largely insensitive to the microturbulence velocity, with individual abundances reduced by only 0.01-0.03 dex going from $\xi_{\rm t}=0$ to $\xi_{\rm t}=1~{\rm km~s}^{-1}$. Microturbulence was not fitted for in the 1D models, as they were used simply for comparative purposes rather than highly accurate determinations. The oxygen and nitrogen abundances used for the 3D, HM and MARCS models were indicative of the most recent values produced with each model (Asplund et al. 2005a, where in the case of oxygen, values are from vibration-rotation and pure rotational OH lines only). The oxygen abundance used for the 1DAV model was derived self-consistently from the same OH lines used to set the oxygen abundances of the other models, and employed previously by Asplund et al. (2004). The nitrogen abundance for the 1DAV model was simply approximated as the 3D value in the absence of any appropriate line calculations. Oxygen and nitrogen abundances used in the line calculations with the different models are shown in Table 1.

Table 1: Oxygen and nitrogen abundances used in line formation calculations for each model in this study. 3D, HM and MARCS values are from Asplund et al. (2004), 1DAV oxygen calculated using the same lines as Asplund et al. (2004) and 1DAV nitrogen abundance simply set to the 3D value.

3 Spectral lines

We utilise six distinct CO line lists, with gf values calculated from Goorvitch & Chackerian (1994). We adopted the 12C16O dissociation energy of Eidelsberg et al. (1987, <)499#>11.108 eV#. All lines were selected to be free of blends. Typical formation heights (estimated for the MARCS model atmosphere) are shown in Fig. 1. CO line formation in higher atmospheric layers is examined in Sect. 5 via bisector analysis of 31 12C16O lines (Table C.1). These are all strong lines that are formed high in the atmosphere. As the height of their formation and hence velocity signatures strongly depend on the temperature structure of the atmosphere, these lines can probe the disputed "COmosphere'' region.

Abundance calculations in Sect. 6 draw upon the five remaining lists, with equivalent widths measured in the ATMOS ATLAS-3 spectrum (cf. Sect. 4). The primary investigation is performed with a set of 13 weak 12C16O lines (Table C.2) of high excitation, as well as sets of 16 13C16O (Table C.3) and 15 12C18O (Table C.4) lines. These lines all have significantly lower formation heights than the strong 12C16O lines considered for Sect. 5 (Fig. 1). Being CO features, in absolute terms they still form reasonably high in the atmosphere, but were selected as forming low enough for abundance determination using the 3D model.

\par\includegraphics[width=7.2cm,clip]{4986f1}\end{figure} Figure 1: Approximate optical depths of core formation for the different CO line lists used in Sect. 5 (strong 12C16O) and Sect. 6 (all others). The temperature structure shown is that of the standard 3D model atmosphere (Asplund et al. 2000b).
Open with DEXTER

The further two lists consist of 15 low excitation (LE) 12C16O (Table C.5) and 66 first overtone ( $\Delta v = 2$) 12C16O (Table C.6) lines. The LE lines are formed across a broad optical depth range high in the atmosphere (due to the low temperatures required for significant populations in the lower energy levels), and the $\Delta v = 2$ at around the same height as the three primary sets (Fig. 1). These supplementary lists are used to derive alternative carbon abundances (and hence isotopic ratios) to the primary (weak) 12C16O list, exploring the abundance performance of the model over higher layers and overtone bands. The three different 12C16O line lists, which also have different temperature sensitivities, have traditionally produced widely varying carbon abundance measures using 1D models, so a consistent result in 3D would greatly increase confidence in the new C abundance and isotopic ratios.

4 The ATMOS infrared solar spectrum

The ATMOS instrument is a Fourier transform spectrograph (FTS). It was carried on a series of space shuttle missions between 1985 and 1994, retrieving pure solar and atmospheric occultation spectra between 625 and 4800 cm-1. Given the instrumental resolution (0.01 cm-1, Farmer 1994) and taking a representative wavenumber of 2000 cm-1 (corresponding to a wavelength of 5000 nm) for the section of the spectrum predominantly considered in this paper, the instrument's resolving power for our purposes is:

\begin{eqnarray*}R = \frac{k}{\Delta k} = \frac{2000}{0.01} = 200~000.

The solar disk-centre spectrum from the 1994 ATMOS ATLAS-3 mission[*] was normalised, with the continuum level defined as the highest intensity point in the extracted section. Given the high signal-to-noise of the data ($\sim$1000:1 around 5000 nm by our measurements), detailed averaging was deemed unnecessary. The Sun's gravitational redshift of 633 m s-1 was removed from the resulting spectrum. The original ATMOS spectrum was apodized using a medium Beer-Norton function (Gunson 2004; Farmer 1994). We hence applied the same apodization to our model output; this essential procedure is described in detail in Appendix A.

To compute CO line bisectors, the absorption profiles were split at their centres, with each side separately interpolated to a 0.01 normalised intensity unit resolution scale using cubic splines. Bisectors were then calculated as the series of wavelength midpoints between the two profile halves. Wavelength errors were estimated conservatively for each bisector point from the previously published (wavelength-dependent) signal-to-noise of the data (Abrams et al. 1996), rather than the higher values we measured from the atlas. Bisectors were truncated above 0.98 normalised intensity, as closer to the continuum they are often dominated by observational noise rather than line asymmetry.

\par\end{figure*} Figure 2: Example spatially and temporally averaged, disk-centre synthesised strong CO line profiles and bisectors (solid lines without error bars), shown in comparison to ATMOS profiles (diamonds) and bisectors (solid lines with error bars). The profile agreement is quite good, although as discussed in the text this requires an unreasonably high C and/or O abundance; the bisectors based on the new 3D model (see text) show better agreement with observation than those of the original 3D model version (dashed lines). The solar gravitational redshift was removed from the ATMOS spectrum, and synthesised profiles have been apodized and fitted in abundance and Doppler shift.
Open with DEXTER

5 CO line shapes and asymmetries

5.1 Analysis and results

In this section, we restrict the discussion to the strongest CO lines, formed in the very high atmospheric layers embroiled in the current temperature debate (e.g. Ayres 2002, and references therein). Given their extreme formation heights, such lines are not ideal abundance determinants, as our study below demonstrates. In Sect. 6, we instead present C and O isotopic abundances based on weaker CO lines formed in slightly lower layers, which our 3D hydrodynamic simulation reproduces more faithfully.

Modelled profiles were fitted in wavelength shift and abundance by minimising the $\chi ^2$ likelihood estimator. The wavelength fitting was permitted to allow for the possibility of systematic errors in the ATMOS wavelength calibration and/or the laboratory wavelengths of the lines considered. Unfortunately, unfitted modelled line strengths were consistently weaker than those observed, requiring that we fit an unrealistically high carbon abundance to produce the best profile agreement. Whilst resultant strong CO line profiles showed reasonable agreement with ATMOS profiles, their bisectors were of a  $\backslash$-shape (Fig. 2, dashed lines) rather than the $\subset$-shape observed (solid lines with error bars).

The hydrodynamic simulations have undergone some slight improvements since the original model atmosphere was generated in 1999, so in an attempt to improve agreement we regenerated the CO line profiles using the most modern version of the simulations available. The newer code contained an improved treatment of the numerical viscosity, and drew on an updated Uppsala package containing slightly improved continuous opacities. The version of the MHD equation-of-state was also more recent, including argon in addition to the 16 other most abundant elements, as well as drawing on slightly altered abundances. The manner in which radiation pressure and energy were included by the MHD tabulation program was also slightly altered. Finally, the mean density and temperature structures used in each case to generate equation-of-state and opacity tables were slightly different, though it is not clear whether this was a primary difference or a secondary effect caused by other updates. The final temperature structures were very nearly identical, however.

Inspection suggested that the failure to properly reproduce strong lines might have been due to incoming gas from the uppermost layer not having time to properly cool and equilibrate to its surroundings before it reached CO line formation heights. To test this hypothesis, a version of the new atmosphere extending to heights of almost 1.2 Mm was also created. The extended simulations were run at a resolution of 50 $\times$ 50 $\times$ 88 and those with the original extension at 50 $\times$ 50 $\times$ 82. The lower resolutions were chosen due to computational time constraints, as the purpose was not to produce the most accurate description of the solar atmosphere available, but to provide a qualitative idea of any improvement in CO bisectors due to the later code or boundary extension. This is an appropriate approach for the qualitative analyses carried out in this Section, but not for e.g. abundance analyses, as carried out in Sect. 6.

Profiles were simulated for all 31 lines, apodized and fitted. Using the regenerated 3D model atmosphere, very good agreement was obtained with observed profiles and bisectors. Unfitted line strengths were still smaller than observed, but the discrepancy was much less than seen previously. Fitted abundances were hence still unrealistically high, though the discrepancy with abundances from weak lines (derived in Sect. 6) was only about half that produced by the original model atmosphere. The scatter in abundances from strong lines was also halved when the regenerated model was employed. Some example resultant profiles and bisectors from the (non-extended) regenerated model atmosphere are shown in comparison to their ATMOS counterparts in Fig. 2. Bisectors derived with the older version of the model atmosphere are also given for comparison. The average reduced $\chi ^2$ value for the agreement between observed and theoretical profiles was an order of magnitude smaller using the later model atmosphere. Profiles and bisectors derived with the extended simulation showed no differences to those produced with the corresponding non-extended version (data not shown). We hence focus mainly on the results from the unextended, later version of the model atmosphere for the remainder of our analysis of line shapes and asymmetries.

5.2 Discussion

To our knowledge, the only prior study of observed asymmetries in solar CO lines is that of Blomme et al. (1994). Without even considering modelled profiles and bisectors, it is plain that the observed ATMOS bisectors do not agree with the $\backslash$-shaped bisectors found previously, a significant result in itself. The reason for this discrepancy appears to simply be that Blomme et al. (1994) extended their bisectors only to intensities of about 0.94 of the continuum intensity, whereas those of Fig. 2 extend to intensities of 0.98. Given that the bisectors of Blomme et al. (1994) were taken from the more noisy SL-3 flight of the ATMOS instrument rather than the ATLAS-3 flight, this lesser extension is not surprising.

The hydrodynamic simulations appear to have been sufficiently improved for some of the deficiencies in the original to be removed. Speculating exactly which minute change or combination of changes described in Sect. 5.1 is responsible for such a higher-order effect as the turning of one end of a bisector is exceedingly difficult. All that can confidently be said on this matter is that it is very unlikely the introduction of argon to the equation-of-state had any effect. Beyond this, any combination of the tiny changes in opacity and abundance inputs, viscosity and radiation pressure treatments or initial mean structure could be responsible.

However, whether the improvement is attributable to resolutional effects rather than actual model differences must be considered. Asplund et al. (2000a) showed that there is an increase in strong Fe line width going from the resolution of the later model to that of the original, with a preferential filling in of the bluewards wing of profiles over the redwards wing. This is due to the asymmetric way in which low resolution models truncate the vertical velocity distribution. This correspondingly manifests itself in a bluewards turning of bisector tops with higher resolution, or equivalently, as a redwards "turning-back'' of bisectors with insufficient resolution.

This "turning-back'' has certainly occurred in the current study, so the important question is not whether, but to what extent the improved bisectors and profiles are attributable to the lower resolution. The bisectors of Asplund et al. (2000a) show a maximum redshifting of the bisector top of 150 m s-1 with the change from a 200 $\times$ 200 $\times$ 82 to 50 $\times$ 50 $\times$ 82 model. They also show a redwards shifting of bisector feet of approximately 100 m s-1, resulting in a net redwards movement of 50 m s-1 of the bisector top relative to the foot. In contrast, the newer model bisectors demonstrate a bare minimum of 250 m s-1 net shift of the top relative to the foot, compared with the (dashed) originals. It would also seem a rather unlikely coincidence for the tenfold difference in $\chi ^2$ values to be due purely to resolutional effects. The improvement of profile shapes with the regenerated model over the original version is also borne out by a much improved convergence in fitted line depths and strengths. Fitted abundances from the 31 lines show far less scatter in the case of the later model, and actual abundance values are also closer to those derived later using weak lines. These changes support the notion that resolutional effects do not play a dominant role in the bisector and profile improvements, as the effects of reduced resolution are expected to be the opposite of that observed, i.e. an increase in scatter (Asplund et al. 2000a) and abundances derived from CO lines (this paper, Sect. 6.4). On the basis of these differences, we estimate the contribution of resolutional effects to be approximately half of the improvement seen.

Whilst agreement has been improved, the final results indicate that the description of the uppermost layers in the model atmosphere is still not perfect. This is not surprising, as many of the physical approximations used begin to break down in the highest layers. Processes such as non-equilibrium hydrogen ionisation, non-LTE continuum formation and extended CO cooling related to non-equilibrium chemistry might actually be required for an accurate description of such regions. Ignoring these non-equilibrium processes could lead to slightly anomalous velocity or temperature structures in the top layers, issues discussed in more depth in Sect. 6.4. Indeed, it turns out that the most likely problem is a slight overestimation of temperature in the uppermost layers of the model, possibly more so in intergranular lanes than granules. This is indicated by the overly high abundances derived from strong and LE lines (as will be seen in Sect. 6.4), as well as our own more extensive internal investigations into the bisector discrepancies and the regeneration exercise (data not shown). Though imperfect, the agreement that was seen with observed bisectors does however indicate that the 3D model at least captures some of the essential physical features of these very high layers, which is encouraging.

5.3 Temperature structure and CO distribution

To use a model to draw accurate conclusions about the temperature structure of the solar atmosphere and the distribution of CO within it, one would ideally require that synthesised CO lines demonstrate excellent agreement with observations from very high atmospheric layers. Unfortunately, despite being the best effort to date, this has not been the case in this study. However, considering that we were able to at least partially reproduce such difficult indicators as the strong line bisectors, and more importantly, that all our investigations point towards the mean temperature structure in the upper layers of the model atmosphere being slightly too high, we can use the existing 3D model atmosphere to indicate temperature upper bounds for the disputed layers. The temperature structure of Fig. 1 hence suggests that cool gas does indeed exist in the lower chromosphere, with an average temperature of <4000 K. Such an upper bound is consistent with the extended 3D LTE atmosphere of Wedemeyer et al. (2004), and the 3500 K proposed by Ayres (2002). Furthermore, the temperature at the site of the previous "minimum'' also seems lower than previously thought, at <4200 K rather than $\sim$4500 K. Owing to the inhomogeneity of the model utilised here, we also tentatively suggest the existence of some persistent gas with T<3700 K at COmospheric heights, and even intermittent gas temperatures of less than 2000 K.

The lack of any profile or bisector improvement in the extended simulations would seem to indicate that the extra layers above 1 Mm probably do not contribute to CO line formation at disc centre for the lines considered here. This also suggests that the downflowing gas entering the upper layer of the non-extended simulations does not play a role in the line formation until it has passed significantly beyond the top boundary layer and been permitted to adjust to the temperature of its surroundings. Hence, gas above a height of $\sim$0.75 Mm might be identified as also generally not contributing to disk-centre solar CO line formation (except perhaps for some extremely strong lines). Furthermore, as it will take some distance to achieve temperature equilibration once gas has passed out of the topmost layer of the simulation, our results could be seen to suggest an uppermost effective extent of the COmosphere of around 700 km, consistent with the indications of Asensio Ramos et al. (2003). It should be noted that this would not however preclude the existence of cool(er) gas at greater heights, just suggests an uppermost extent of CO in significant concentrations.

6 C and O isotopic abundance determinations

6.1 Dealing with CO and isotopes

Working with a known oxygen abundance, a fitted abundance difference for CO lines can be approximately interpreted as a difference in solar carbon abundance, even though strictly speaking it is actually indicative of a difference in CO abundance at the height of line formation. In order to use CO lines to derive a solar carbon abundance, the input carbon abundance was therefore iteratively altered according to the average fitted abundances of the CO lines in the previous iteration, until convergence was obtained. Oxygen concentration effects the equilibrium position of the CO formation-dissociation reaction, so carbon abundances derived in this manner are dependent upon the adopted oxygen abundance. Input oxygen abundances were hence carefully chosen (cf. Sect. 2) and kept constant throughout the study.

Using the above technique, the solar 12C abundance was determined separately using three different 12C16O line lists. Since the 3D model input abundances do not differentiate between isotopes, these determinations were performed by altering the total input carbon abundance, though the final abundance arrived at is not overall abundance, but the carbon abundance if all C in CO were contained in 12C16O. This can approximately be called the 12C abundance, even though it does not include any 12C tied up in 12C17O or 12C18O. Though they are minimal given the high values of 16O/18O and 16O/17O, these contributions were taken into account in the calculation of a total carbon abundance by the use of a fractional scalefactor (Appendix B.3).

The concentrations of 13C16O and 12C18O relative to 12C16O are identical to the 12C/13C and 16O/18O ratios. In order to determine the concentrations of the CO isotopomers, given the lack of provision for isotopic differentiation in the model, the isotopomeric lines were treated in the radiative transfer simulations as if they were created by 12C16O. In order to do this correctly, the radiative transfer code had to be altered to include provision for a mass scalefactor and an opacity scalefactor (refer to Appendix B).

6.2 Abundance calculations

For each iteration, three profiles with $\log gf$ values differing by 0.2 dex were calculated for a given line and interpolated between using cubic splines to arrive at the desired line strength. As described in Sect. 2, abundance calculations were carried out for all five line lists with four different models: 3D, HM, MARCS and 1DAV. The 3D model used was that of Asplund et al. (2000b), the same model atmosphere used for all previous 3D abundance determinations (e.g. Asplund 2000; Asplund et al. 2005b; Asplund 2004; Asplund et al. 2004,2000c,2005a). Since accuracy is paramount in abundance determinations, the regenerated version described in Sect. 5 would have been inappropriate given its reduced resolution. This model was also used over the regenerated version for consistency and comparability with the previous abundance determinations.

Both equivalent width and profile fitting (the latter via $\chi ^2$-analysis) were tested upon the first iteration of the weak 12C16O lines, with virtually no difference found in calculated abundance. However, some of the lines could not be effectively profile fitted by minimising the $\chi ^2$ statistic without some form of masking, due to the presence of other nearby lines in the ATMOS spectrum. For this reason, equivalent width fitting was used for all subsequent abundance measures, as this problem was only likely to worsen when the weaker isotopomeric lines were considered. The same Beer-Norton apodization as used in Sect. 5 (medium BNA, characteristic velocity 1.5 km s-1) was also applied to each line, though more for the sake of consistency than anything else. Convolution with any function normalised to unit area (as the medium Beer-Norton used was) should be an area-preserving operation, thereby not effecting line equivalent widths nor therefore the abundances determined with them. For completeness however, it should be noted that in the case of the $\Delta v = 2$ lines, a BNA characteristic velocity of 0.75 km s-1 was used, reflecting the higher resolving power of the ATMOS instrument at the shorter wavelengths of these lines (see Table C.6).

\mbox{\includegraphics[width=8cm]{4986f3a}\hspace{5mm} \includeg...
...{4986f3b}\hspace{5mm} \includegraphics[width=8.2cm]{4986f3d} }
\par\end{figure*} Figure 3: Solar carbon abundances indicated by the weak 12C16O lines, displayed according to equivalent width ( top) and excitation potential ( bottom). On the left, filled circles indicate 3D results and open circles 1DAV results, whilst on the right filled circles are HM values and open circles MARCS results. Trendlines are produced as linear fits to data sets using a minimised $\chi ^2$ method placing equal weight on each point, with solid lines corresponding to filled circles and dashed lines to open circles. No significant trends can be seen in the output of any model.
Open with DEXTER

...4986f4b}\hspace{5mm}\includegraphics[height=52.5mm]{4986f4d} }
\par\end{figure*} Figure 4: The same as Fig. 3, but using the LE 12C16O lines. Definite trends can be seen with equivalent width in the output of the 3D and MARCS models. Significantly, the 12C abundances implied by the weakest (i.e. lowest formation height) lines in 3D are consistent with the abundances derived using the weak and $\Delta v = 2$ 12C16O lines. The agreement deteriorates with increased LE 12C16O line strength and therefore formation height. Less prominent trends are also evident in excitation potential for all models.
Open with DEXTER

...4986f5b}\hspace{5mm}\includegraphics[height=52.5mm]{4986f5d} }
\par\end{figure*} Figure 5: The same as Fig. 3, but using the $\Delta v = 2$ 12C16O lines. Significant trends can be seen with excitation potential in the 1D results, though none in the case of the 3D model.
Open with DEXTER

...4986f6b}\hspace{5mm}\includegraphics[height=52.5mm]{4986f6d} }
\par\end{figure*} Figure 6: The same as Fig. 3, but indicating 13C through the use of the 13C16O lines. No significant trends can be seen in the output of any model.
Open with DEXTER

...4986f7b}\hspace{5mm}\includegraphics[height=52.5mm]{4986f7d} }
\par\end{figure*} Figure 7: The same as Fig. 3, but indicating 18O through the use of the 12C18O lines. Significant trends can be seen in equivalent width and excitation potential in the HM and MARCS results.
Open with DEXTER

6.3 Results

Derived abundances are plotted as a function of equivalent width and excitation potential for each line in the weak 12C16O (Fig. 3), LE 12C16O (Fig. 4), $\Delta v = 2$ 12C16O (Fig. 5), 13C16O (Fig. 6) and 12C18O sets (Fig. 7). The only appreciable trend evident in the 3D results is with equivalent width in the case of the LE lines, though these lines do also produce a smaller trend in excitation potential in 3D. The HM and MARCS models display trends in equivalent width and excitation potential for the 12C18O lines, and excitation potential in the $\Delta v = 2$ lists. In addition, the MARCS model exhibits a trend with equivalent width in its LE results. In the case of the 1DAV model, no significant trends are seen except for a small slope with equivalent width in the 12C18O list.

Final abundance measures and isotopic ratios as defined by the weak, LE and $\Delta v = 2$ 12C16O lines are tabulated (Table 2). Errors in 12C abundances are given by single standard deviations within individual 12C16O line lists. Errors in bulk C abundances are single standard deviations within each 12C16O list when a fractional scalefactor (cf. Appendix B.3) was employed, as uncertainty in the scalefactor due to isotopic abundance uncertainties was negligible compared to line-to-line scatter. Similarly, errors in isotopic abundances are a single standard deviation within the relevant isotopomeric list. Errors in isotopic ratios were determined in logarithmic units as the sum in quadrature of errors attached to logarithmic 12C and isotopic abundances. Asymmetric errors in absolute ratios are given as corresponding deviations due to these calculated (symmetric) logarithmic errors. Overall results in 3D from different indicators are summarised in Table 3.

6.4 The solar C abundance

The absence of significant trends in 3D model abundances over the primary (weak) and overtone ( $\Delta v = 2$) line lists (Figs. 3 and 5), as well as the similarity of the resulting average abundances ( $\log\epsilon_{\rm C}$ = 8.40, 8.37), give us confidence in the accuracy of these values. They also suggest that the 3D model is accurate in the quite high atmospheric layers in which these lines form. The agreement of abundances derived using the normal weak and overtone lines is quite an achievement, as these indicators have not generally produced consistent results. The HM model for example produces abundances differing by 0.08 dex between the two lists (and by nearly three times this between CO lines and other indicators of carbon abundance stated in Asplund et al. 2005b). It also exhibits a clear trend with excitation potential for the $\Delta v = 2$ lines (bottom-right of Fig. 5), further suggesting the inadequacy of its description of overtone CO. The same trends are true of MARCS, though it at least produces consistent abundances between the two line lists. That the 1DAV model abundances in these two cases are much closer to the 3D results, though still not quite as low, reflects the fact that both the mean temperature structure and the temperature inhomogeneities play a role in producing different CO-derived C abundances in 3D than with HM or MARCS.

The same success was not seen in LE results, where a striking trend with equivalent width is present in the 3D case (top left of Fig. 4). The reason for this is very likely the larger range of formation depths of these lines and increasing inadequacies of the 3D model atmosphere in higher layers. The weaker of the LE lines can be seen to produce abundance measures very close to the primary and overtone results, as they form at similar heights (cf. Fig. 1). However, the derived abundance increases with greater equivalent width and line formation height, resulting in an average abundance quite a bit higher than the primary or overtone diagnostics.

Table 2: Abundances and isotopic ratios implied by the weak ( $\Delta v = 1$), LE ( $\Delta v = 1$) and $\Delta v = 2$ 12C16O line lists as indicator of 12C abundance. Note the large differences between 1D and 3D isotopic ratios. Note also that the HM and MARCS 1D models indicate higher abundances than the 1DAV, which in turn produces a higher abundance than the 3D model.

Table 3: Summary of carbon abundances and isotopic ratios produced with the 3D model by the three different 12C16O line lists, as well as the adopted values. Note the poor agreement of the LE results with the other two, reflecting the trend seen in Fig. 4. Adopted values were calculated via 2:0:1 weightings of the weak:LE: $\Delta v = 2$ lists.

A number of reasons could be postulated for this discrepancy. The failure of the chemical equilibrium approximation (ICE) with height in the atmosphere would mean that CO density was being overpredicted in the models, producing stronger lines than otherwise would result from a particular abundance. Hence, less abundance would be required to reproduce a given line profile than in the non-ICE case, not more as is seen here. The breakdown of LTE at height in the atmosphere could possibly cause overestimated abundances here, though it seems unlikely given the repeated conclusion that CO lines form in LTE (Ayres & Wiedemann 1989; Uitenbroek 2000b). Another possibility could be that temperature contrast in the upper layers of the model is too low, as an increase in this contrast would produce lower temperature cool regions, which contribute more to increasing line strength than hotter regions would to decrease it, due to the increase in CO and the nonlinear temperature dependence of line formation. However, this seems unlikely as our bisector investigations hinted at too much temperature contrast high in the atmosphere.

Also possible is that the velocity structure of the upper atmosphere is slightly incorrect in the 3D simulations. The MARCS model exhibited a similar equivalent width trend in implied abundances to the 3D model for the LE lines, also suggesting problems in its upper layers. However, even the strongest LE lines were not sensitive to the microturbulence parameter used with the MARCS model, making this explanation unlikely. Similarly, other controls revealed that the problem was not due to the adopted collisional damping parameter. The most likely explanation for the pronounced trend seen in the 3D LE results is that the mean temperature structure at height is slightly too high, causing a reduction in line strength for a given abundance and conversely, increased abundances for given line strengths. This is consistent with the lack of trend but higher average abundance indicated by the 1DAV model, as the horizontal averaging would smear out the effects of overly hot localised regions (e.g. intergranular lanes) at height but produce higher abundances overall due to the complete absence of temperature inhomogeneities. This is also in line with the bisector discrepancies seen using strong CO lines.

The regenerated model of Sect. 5 was trialled on the LE and weak 12C16O lines to see if agreement could be improved, and also to check that the newer model did not produce significantly different abundances where the older model was thought to be accurate (i.e. the weak lines). Derived abundances in both cases were approximately 0.03 dex larger, virtually identical to the effect of low resolution upon iron abundances found by Asplund et al. (2000a).

Despite the poor performance of the 3D model for the LE lines, none of the 1D models managed to derive much better agreement between these and the weak 12C16O lines (which are regarded as the best CO indicators of carbon abundance). In the sense of the LE lines' disagreement with other diagnostics, the 3D model simply presented no improvement, rather than any loss in performance over HM and MARCS. In all three diagnostics, the carbon abundance in 3D is considerably lower than indicated by HM or MARCS, as has generally been found for other species (e.g. Asplund et al. 2005b,2004,2000c,2005a). This is a general consequence of extending atmospheric simulations to three dimensions, due to the permission of temperature inhomogeneities and the nonlinear temperature dependence of line formation. That all three sets of 12C16O lines continued to display this effect is a positive comment on the accuracy and consistency of the 3D model, especially given the past difficulties with 1D analyses of CO lines (Grevesse et al. 1995).

Given the relative performance of the different 12C16O line lists and their previously recognised suitability for carbon abundance determination, for the final analysis weak lines were given a weighting of 2, $\Delta v = 2$ lines a weighting of 1 and LE lines no weighting at all. We note however that abundances from the weakest LE lines are in good agreement with those from the two preferred line types. Hence, on the basis of the lines measured in this study, the final carbon abundance arrived at is

\begin{eqnarray*}\log\epsilon_{\rm C} = 8.39\pm0.05.

The stated error is designed to encompass unquantified systematic errors from e.g. atomic/molecular data, the model atmosphere and equivalent width measurements; these systematics almost certainly outweigh any statistical variability amongst abundances indicated by individual lines or lists.

6.5 The solar 12C / 13C and 16O / 18O ratios

The abundances produced by the individual 13C16O and 12C18O lines in 3D (left panels in Figs. 6 and 7) show no significant equivalent width or excitation potential dependence, hence indicating no major inadequacies in the relevant atmospheric layers of the model. This is also the case for 1D calculations with the 13C16O lines (right panels in Fig. 6) but not the 12C18O lines, where trends with equivalent width and excitation potential are evident in the HM and MARCS results (right panels in Fig. 7). The two trends are explicitly linked, as the higher excitation lines form lower and therefore have smaller equivalent widths, so it is no surprise that one trend exists in equivalent width and the opposite is seen in excitation potential. These trends would seem to indicate the superiority of the 3D model in the description of the very weak 12C18O lines.

The derived isotopic ratios in Table 2 using the different model atmospheres show a very strong separation between 3D and 1D results. This is surprising at first, though suggests that in this context, the move to three dimensions is an important improvement. This is highlighted in particular by the general difference between 3D and 1DAV results. Our investigations have revealed that, due their lower excitation potentials, the isotopomeric lines are actually over 20% more temperature sensitive than the weak 12C16O lines. At least in the case of the weak 12C16O lines, this is almost certainly the reason for the greater decrease in isotopic abundances than the carbon abundance in 3D. This is because the heterogeneic temperature structure of the 3D simulations (which is responsible for reduced abundances normally) has a greater decreasing effect upon the isotopomeric lines and therefore 13C and 18O abundances than it does upon the derived carbon abundance. Because of this differential temperature sensitivity, the horizontal averaging implicit in 1D models means that they have difficulty reproducing the correct solar isotopic ratios. Given their excitation potentials, the temperature sensitivities of the LE and $\Delta v = 2$ lines are rather similar to the isotopomeric lines. This is the reason for the better agreement between 1D and 3D isotopic ratios from the $\Delta v = 2$ lines than from the weak 12C16O lines. Any similar effect in ratios derived from LE lines is lost in the distortion of the 3D carbon abundance by the observed trend with formation height.

Because each 12C16O list produced different abundances and each was compared to the same 13C16O and 12C18O results, different 3D isotopic ratios (and absolute 18O abundances, since $\epsilon_{\rm O}$ was fixed) were obtained. However, the LE ratios can be discarded due to the already established invalidity of the LE carbon abundance. In the derivation of the final isotopic ratios, the $\Delta v = 2$ ratios were given a lower weighting in the same manner as in the calculation of the final carbon abundance (i.e. half the weighting afforded the weak ratios). Hence, the adopted isotopic ratios are:

\begin{eqnarray*}&& ^{12}{\rm C}/^{13}{\rm C} = 86.8^{+3.9}_{-3.7}\\
&& ^{16}{\rm O}/^{18}{\rm O} = 479^{+29}_{-28}.

Errors are combinations of statistical errors in the weak and $\Delta v = 2$ results. We assume that the systematic errors discussed as relevant to the carbon abundance in Sect. 6.4 effect all lines considered for ratios equally, so need not be included in final ratio errors. This assumption is based on fact that the isotopic, weak and $\Delta v = 2$ lines all form at approximately the same photospheric heights, reflecting essentially the same line formation processes.

7 Comparisons with previous work

This study represents a significant step forward from the only other 3D investigation of CO line formation to date (Uitenbroek 2000a). Firstly, we used a more modern 3D hydrodynamic model atmosphere code, whereas the previous study was based upon the earlier Stein & Nordlund (1989) version, running at a resolution of just 63 $\times$ 63 $\times$ 63. Secondly, the current study produced profiles spatially averaged over the entire simulation domain and temporally averaged over about 100 snapshots corresponding to approximately 50 min of solar time, whereas that of Uitenbroek (2000a) used a vertical slice through a single snapshot. Thirdly, the current study extended about 1 Mm above the solar surface, compared with 0.6 Mm in the earlier work. Finally, agreement between the modelled spectrum and observation was excellent overall, with reasonable agreement even found at the level of individual line bisectors, a measure sensitive to small deviations of the simulation from reality. As suggested by Asensio Ramos et al. (2003), the overly cool version of the 3D model atmosphere employed by Uitenbroek is probably somewhat responsible for the overly deep CO line cores he found. The recently reduced solar carbon and oxygen abundances (Asplund et al. 2005b,2004) have also played a role in the success of the current study, permitting the deviation of line profiles from the Grevesse & Sauval (1998) abundances in order to improve agreement with observation. This option was not realistically available to Uitenbroek (2000a) given prevailing wisdom at the time, and in light of the current results this probably also contributed to his modelled overprediction of CO line depths.

In the context of spectral line formation, our results also constitute an improvement over the more recent 2D work of Wedemeyer-Böhm et al. (2005). This is due not only to the better agreement achieved with observation, but also the superiority of our model atmosphere in this context. Apart from being 3D rather than 2D and therefore inherently more realistic, the model we employ includes line blanketing, rather than treating radiative transfer in the hydrodynamic code as grey. In addition, our model exhibits an effective temperature in excellent agreement with the Sun's (see Asplund et al. 2000b), whereas the snapshots Wedemeyer-Böhm et al. employ have an effective temperature approximately 380 K lower. Whilst Wedemeyer-Böhm et al. utilise non-equilibrium chemistry, they find that it is not required at the heights of formation considered in this paper, justifying our use of the ICE approximation throughout.

The adopted bulk carbon abundance is in excellent agreement with the $\log\epsilon_{\rm C}=8.39$ $\pm$ 0.05 of Asplund et al. (2005b), which was based on very different indicators (C I, [C I], C2 and CH). The current figure therefore firms our belief in the accuracy of both carbon results and the 3D model atmosphere. In comparison, CO results with the preferred lines using the HM and MARCS models do not agree at all with their C I, [C I], C2 and CH counterparts given by Asplund et al. (2005b); CO lines give $\log\epsilon_{\rm C} = 8.60$, 8.69 but other lines 8.39-8.53 with HM, whilst with MARCS CO lines indicate $\log\epsilon_{\rm C} = 8.55$, 8.58 and other lines give 8.35-8.46. Given the height of formation of the CO lines, this result is a remarkable success for the 3D model. The new carbon abundance constitutes a major reduction from the commonly adopted Grevesse & Sauval (1998) value of $\log\epsilon_{\rm C} = 8.52$ $\pm$ 0.06, not to mention a multitude of even higher earlier estimates (e.g. Harris et al. 1987; Anders & Grevesse 1989).

Using the HM model, Harris et al. (1987) found isotopic ratios of 12C/13C = 84 $\pm$ 5 and 16O/18O = 440 $\pm$ 50, broadly consistent with the ratios we obtain according to the stated errors. In comparison, our own HM results (weighted amongst the lists in the same manner as the 3D ratios) were 12C/13C = 73.4 +2.8-2.7 and 16O/18O = 368 +25-23. The difference in HM results between the two studies evidently reflects some combination of the improved observations, line lists and molecular data used in the present work[*]. However, the effect of these improvements is to reduce the ratios, whereas the adopted ratios calculated in 3D are slightly greater than those previously found. Hence, the agreement between our results and those of Harris et al. (1987) is somewhat fortuitous; any difference between old and new adopted ratios primarily reflects the shift from 1D to 3D modelling, with the effects of this shift in fact mostly nullified by improved input data.

We also believe that our ratios are more reliable than those very recently produced independently by Ayres et al. (2006; 12C/13C = 80 $\pm$ 3 and 16O/18O = 440 $\pm$ 20 using their preferred model, where errors are standard deviations rather than the lower standard error Ayres et al. favour). They utilise our 1DAV model and a multicomponent version of it meant to approximate the true 3D case. Discrepancies between isotopic ratios derived in the two studies using the 1DAV model can be explained by differences in measured equivalent widths and our use of the opacity scalefactors detailed in Sect. B.2[*]. We believe our equivalent width measurements to be superior: whilst Ayres et al. (2006) utilise Gaussian fits to ATMOS data, we performed direct integration upon line profiles, since the ATMOS ATLAS-3 data is of sufficiently high signal-to-noise to make this possible and observed profiles always exhibit more extended wings than a pure Gaussian. Furthermore, all Ayres et al.'s Gaussian fits to isotopomeric lines are constrained to have a FWHM of 4.05 km s-1, whereas our integration allowed this value to be directly determined. However, Ayres et al. choose to discard the results of our model in favour of an empirical 1D model based upon a "Double-Dip'' photospheric structure, which has no corroborating physical or theoretical basis beyond its ability to match the observed centre-to-limb variation of CO lines and the solar continuum. The derived abundances show marked trends in excitation potential and a disagreement between fundamental and overtone bands, which our 3D and DAV models do not (either in this paper or that of Ayres et al. 2006). To compensate for these trends, the authors introduce a highly-questionable, excitation potential-dependent change in given gf values. Ayres et al.'s choice of model is made upon the basis of the 1DAV model failing to properly reproduce the observed continuum and CO core centre-to-limb behaviour. However, the stated discrepancy in the continuum centre-to-limb results is about double that produced by our own calculations. We intend to return to this important issue in a subsequent study, where the temperature structure and full results of a large number of observational tests will be described.

The 12C/13C ratio we find is in excellent agreement with the Coplen et al. (2002) terrestrial value of 89.4 $\pm$ 0.2. Our ratio corresponds to a  $\delta^{13}$C value of 30+46-44 (without taking into account the uncertainty in the $\delta$-scale), in agreement with some of the lunar regolith measures discussed by Wiens et al. (2004). Our ratio does not agree with the $\delta^{13}{\rm C}\leq -105$ $\pm$ 20 measured by Hashizume et al. (2004), suggesting that theories for the formation of organics in the early solar system based upon this low value may require revision. We confirm a high solar 12C/13C ratio relative to the local ISM (62 $\pm$ 4: Langer & Penzias 1993; 68 $\pm$ 15: Milam et al. 2005). This may support the notion of a galactic enrichment of 13C over the Sun's lifetime (e.g. Milam et al. 2005). Alternatively, the difference could be due to greater isotopic inhomogeneity in the galaxy than previously thought, and migration of the solar system through it.

Our derived 16O/18O abundance is also consistent with the terrestrial ratio of 498.7 $\pm$ 0.1 (Coplen et al. 2002). The corresponding $\delta^{18}$O value is 41+67-59 (with error in the $\delta$-scale not accounted for). Despite the lack of error bars on the $\delta^{18}{\rm O} = -50$ prediction of the self-shielding model (Yurimoto & Kuramoto 2004; Yin 2004), it seems unlikely that our results would be consistent with this prediction. However, our result is in excellent agreement with the $\delta^{18}{\rm O} = 50$ recently determined by Ireland et al. (2006) from lunar grains irradiated by the solar wind. However, given uncertainty about how well the solar wind actually reflects the solar composition, and the size of our estimated errors, this amazing agreement could be slightly fortuitous.

8 Conclusions

This is the most successful modelling of CO line formation in the solar atmosphere to date. Solar abundance determinations using CO lines produced a carbon abundance of $\log\epsilon_{\rm C}=8.39$ $\pm$ 0.05, and isotopic ratios of 12C/13C = 86.8 +3.9-3.7 and 16O/18O = 479 +29-28. These results represent a significant improvement over those of the past due to the combination of a state-of-the-art convective 3D model atmosphere, updated atomic data, better line lists and more accurate observations. The carbon abundance is in excellent agreement with the recent findings based upon entirely different indicators of Asplund et al. (2005b), suggesting that the past problems with CO-derived abundances have been solved. Both ratios are in excellent agreement with corresponding telluric values, though the oxygen ratio is in even closer agreement with a non-terrestrial value inferred by the latest lunar regolith analysis (Ireland et al. 2006). Our results also support the existence of a COmosphere, with a representative temperature below 4000 K.

We dedicate this article to the memory of Hartmut Holweger, who passed away in early April, 2006. He made many important contributions to solar (and stellar) physics, including studies of the solar chemical composition and the atmospheric temperature structure, both topics of this work. He will be deeply missed by his many colleagues and friends.
We are grateful to Regner Trampedach for assistance with computing 3D solar models and useful discussions, Mike Gunson for ATMOS advice and Charles Jenkins for statistical assistance. Thanks to Tom Ayres for helpful discussions and providing the preprint of his latest paper, and to Frank Briggs and Nathan Kilah for further discussions. P.S. acknowledges support from the ANU National Undergraduate Scholarship, M.A. from the Australian Research Council and N.G. from the Royal Belgian Observatory. This work has made extensive use of the Astrophysical Data System.



Online Material

Appendix A: ATMOS apodization

Knowing the instrument's resolving power, it becomes possible to simulate the observational effect of an FTS. This is simply done by convolving modelled spectra with an instrumental line shape (ILS) characterised by the spectral resolution in the appropriate scale.

Unfortunately however, FTS instruments are not always so well behaved as to exhibit perfect top-hat interferogram profiles. This was certainly the case with the ATMOS instrument, which had an ILS that was far from a perfect sinc (Gunson 2004). To deal with the uncertainty in the true instrumental profiles of transform spectrometers, apodization is sometimes employed. The process of apodization involves the multiplication of the interferogram with a profile other than the typical top-hat function, such as a Gaussian or simple triangle function, with somewhat tapered edges. This of course has the same effect as convolving the spectral output with an ILS that is the Fourier transform of the apodizing function which one multiplies the interferogram by. This smears out the variability in the instrumental profile, and dampens harmonic edge effects produced in the resultant spectrum by the finite maximum optical path difference of the FTS.

The apodization applied by the ATMOS team (Gunson 2004; Farmer 1994) was the medium function of Norton & Beer (1976):

                              $\displaystyle %
ILS(\sigma)$ = $\displaystyle 0.26{\rm sinc}~a - 0.464514\frac{{\rm sinc}~a - \cos
    $\displaystyle - 13.422570\frac{(1-3/a^2){\rm sinc}~a + (3/a^2)\cos a}{a^2},$ (A.1)

where $a = \frac{\pi\sigma}{\Delta\sigma}$, $\sigma$ is the spectral variable and  $\Delta \sigma $ its resolution. This function is one of three purposefully created by Norton & Beer for their minimal resolutional broadening of spectral points and maximal damping of secondary maxima. In order to correctly simulate the effects of the ATMOS instrument and the post-processing performed by its operators, we convolved modelled profiles with the same medium Beer-Norton function. Model spectra were expressed in terms of Doppler velocity, so the convolving function argument was velocity resolution, i.e. $\frac{c}{R}$ = 1.5 km s-1. The striking effect of this post-processing is seen in Fig. A.1, where a predicted 3D line profile that has undergone Beer-Norton apodization (BNA) fits ATMOS data far better than either its unconvolved counterpart, or a profile convolved with a plain sinc function.

\par\includegraphics[width=7.6cm,clip]{4986fA1}\end{figure} Figure A.1: The effects of medium BNA upon the same 4752.2 nm CO line as shown in Fig. 2. The dotted profile has had no instrumental profile applied to it, whereas the solid line illustrates the effects of medium BNA characterised by $\Delta \sigma $ = 1.5 km s-1. Convolution with a simple sinc function of $\Delta \sigma $ = 1.5 km s-1 leaves the line profile virtually identical to its unconvolved counterpart. In the interests of clarity this curve is not shown. Note the far better agreement with the ATMOS spectrum of the apodized than the unapodized profile (reduced $\chi ^2$ measures indicate an order of magnitude difference), vividly demonstrating the importance of correctly emulating instrumental effects when working with this data.

Appendix B: Scalefactors

B.1 Mass scalefactor

The mass scalefactor allowed the correct molecular mass to be used in the line formation calculations, and was simply the ratio of the mass of the isotopomer in question to that of the most common isotopomer (i.e. 12C16O). These factors, 1.03583 for 13C16O and 1.07157 for 12C18O, were calculated from nuclear weights given by Audi & Wapstra (1995).

B.2 Opacity scalefactor

The opacity scalefactor contained the information about different isotopomer concentrations and therefore the atomic isotopic ratios, as the line opacities for any particular transition in two different isotopomers are proportional to their densities. The use of an opacity scalefactor was preferable to simply altering the input carbon abundance in order to emulate the difference in isotopic abundances. This is because the overall carbon abundance was not altered and therefore did not feed back on other parts of the simulation like CN line formation or the temperature structure, which might in turn indirectly effect CO line formation. The opacity scalefactor was necessary anyway, in order to include an actual molecular difference (rather than just that due to abundances) in opacity between isotopomers. This opacity correction arises as follows:

First, recall that Doppler (thermal) line broadening by a species of mass m has a narrow Gaussian characteristic profile of the form

\phi_{\rm D}(\Delta\lambda) = \frac{1}{\sqrt{\pi}\Delta\lambda_{\rm D}}{\rm e}^{-(\Delta\lambda/\Delta\lambda_{\rm D})^2},
\end{displaymath} (B.1)

where $\Delta\lambda$ is the distance from the line centre of wavelength $\lambda_0$ and

\Delta\lambda_{\rm D} = \frac{\lambda_0}{c} \sqrt{\frac{2kT}{m} + \xi_{\rm t}^2}
\end{displaymath} (B.2)

is the Doppler width of the Gaussian. The microturbulent velocity  $\xi_{\rm t}$ is an extra "fudge factor'' used in analyses based on 1D model atmospheres, introduced to emulate Doppler broadening due to inherently three-dimensional microturbulence in the gas.

Now, consider two isotopomers of a diatomic molecule of elements X and Y: " $\mathcal{A}$'' = aXcY and " $\mathcal{B}$'' = bXcY. The opacity of a species in a given transition depends on the density of absorbers, the transition probability and the relevant line broadening effects. In the case of a weak lines, as those considered in this study are, Doppler broadening given by the profile  $\phi_{\rm D}$ from Eq. (B.1) can be approximated as the only significant contributor to the latter, meaning that for our molecule  $\mathcal{B}$,

\kappa_\mathcal{B^\star} \varpropto N(\mathcal{B}^\star)f_{\mathcal{B}^\star}\phi_{\rm D}
\end{displaymath} (B.3)

where $\mathcal{B}^\star$ denotes the isotopomer  $\mathcal{B}$ in some given energy level from which the transition in question occurs, $N(\mathcal{B}^\star)$ refers to the number density of isotopomers  $\mathcal{B}$ in the excited energy level and  $f_{\mathcal{B}^\star}$ refers to the oscillator strength or probability of the particular transition from the excited energy level. Now, at any given wavelength width  $\Delta\lambda$ in Eq. (B.1), we see that
    $\displaystyle \phi_{\rm D} \varpropto \frac{1}{\Delta\lambda_{\rm D}}\hspace{2pt}\varpropto\hspace{2pt}\frac{1}{\sqrt{\frac{2kT}{m} + \xi_{\rm t}^2}}$  
    % latex2html id marker 8157
$\displaystyle \therefore \kappa_\mathcal{B^\star} \...
...r)f_{\mathcal{B}^\star}}{\sqrt{\frac{2kT}{m_\mathcal{B}} + \xi_{\rm t}^2}}\cdot$ (B.4)

Now, assuming LTE, we know from statistical mechanics (Mihalas 1978) that a Boltzmann distribution of energy levels implies that where many instances of some microscopic system $\Upsilon$ occur, the population of a particular excited energy level  $\Upsilon^\star$ is given by

N(\Upsilon^\star) = \frac{N(\Upsilon)g_{\Upsilon^\star}\exp(\frac{-E_{\Upsilon^\star}}{kT})}{Z(\Upsilon)},
\end{displaymath} (B.5)

where $g_{\Upsilon^\star}$ is the statistical weight of the excited energy level, $E_{\Upsilon^\star}$ is its energy and  $Z(\Upsilon)$ is the partition function of the system. Hence,

\kappa_\mathcal{B^\star} \varpropto
...hcal{B})\sqrt{\frac{2kT}{m_\mathcal{B}} + \xi_{\rm t}^2}}\cdot
\end{displaymath} (B.6)

Now, the number of molecules is related to the number density of the molecule's constituent atoms through the well-known Guldberg-Waage law, so

N(\mathcal{B}) = \frac{N(^b{\rm X})N(^c{\rm Y})}{\varkappa(\mathcal{B})},
\end{displaymath} (B.7)

where $\varkappa(\mathcal{B})$ is the equilibrium constant of the formation-dissociation reaction of the isotopomer  $\mathcal{B}$, given (Schadee 1964; Tatum 1966) by

\varkappa(\mathcal{B}) \varpropto \frac{Z(^b{\rm X})Z(^c{\r...
\end{displaymath} (B.8)

where $D_\mathcal{B}$ and $\mu_\mathcal{B}$ are the dissociation energy and reduced mass respectively of isotopomer  $\mathcal{B}$. Hence, combining Eqs. (B.6)-(B.8) we now see that
                                       $\displaystyle \kappa_\mathcal{B^\star} \varpropto
\frac{N(^b{\rm X})N(^c{\rm Y}...
...-E_{\mathcal{B}^\star}}{kT})}{\sqrt{\frac{2kT}{m_\mathcal{B}} + \xi_{\rm t}^2}}$  
    $\displaystyle \hspace{6mm} \varpropto
\frac{N(^b {\rm X})N(^c {\rm Y})Z(\mathca...
...\rm X})Z(^c {\rm Y})\mu_\mathcal{B}^\frac{3}{2}\exp(\frac{-D_\mathcal{B}}{kT})}$  
    $\displaystyle \hspace{9mm}\times
\frac{gf_{\mathcal{B}^\star}\exp(\frac{-E_{\mathcal{B}^\star}}{kT})}{\sqrt{\frac{2kT}{m_\mathcal{B}} + \xi_{\rm t}^2}}$  
    % latex2html id marker 8176
$\displaystyle \therefore\hspace{1cm}\kappa_\mathcal...
...})\mu_\mathcal{B}^\frac{3}{2}\sqrt{\frac{2kT}{m_\mathcal{B}} + \xi_{\rm t}^2}},$ (B.9)

where $g_{\mathcal{B}^\star}$ and $f_{\mathcal{B}^\star}$ have been combined into  $gf_{\mathcal{B}^\star}$, the gf value of the line.

So, the opacity scalefactor $\varsigma_\kappa$ for an isotopomer $\mathcal{B}$, where the reference isotopomer is $\mathcal{A}$ is then

                        $\displaystyle %
\varsigma_\kappa \equiv \frac{\kappa_\mathcal{B^\star}}{\kappa_\mathcal{A^\star}}$ = $\displaystyle \frac{\frac{N(^b{\rm X})N(^c{\rm Y})gf_{\mathcal{B}^\star}\exp(\f...
...})\mu_\mathcal{A}^\frac{3}{2}\sqrt{\frac{2kT}{m_\mathcal{A}} + \xi_{\rm t}^2}}}$  
  = $\displaystyle \frac{N(^b{\rm X})gf_{\mathcal{B}^\star}Z(^a{\rm X})\mu_\mathcal{...
    $\displaystyle \times
\frac{\sqrt{\frac{2kT}{m_\mathcal{A}} + \xi_{\rm t}^2}}{\sqrt{\frac{2kT}{m_\mathcal{B}} + \xi_{\rm t}^2}}\cdot$ (B.10)

Now, we make the approximations $gf_{\mathcal{B}^\star} \approx gf_{\mathcal{A}^\star}$, $D_\mathcal{B} \approx D_\mathcal{A}$ and $E_{\mathcal{B}^\star} \approx E_{\mathcal{A}^\star}$, since  $\mathcal{A}$ and  $\mathcal{B}$ are very nearly the same molecule so their transitions and energy levels will be almost identical. Further, assuming $Z(^b{\rm X}) \approx Z(^a{\rm X})$ due to the first-order dependence of atomic orbitals upon charge (and not mass), we arrive at

\varsigma_\kappa = \frac{N(^b{\rm X})\mu_\mathcal{A}^\frac{...
...^\frac{3}{2}\sqrt{\frac{2kT}{m_\mathcal{B}} + \xi_{\rm t}^2}},
\end{displaymath} (B.11)

where the final scalefactor is independent of transition. The reduced mass opacity corrections used for 13C16O and 12C18O were therefore 0.9346 and 0.9293 respectively, where nuclear masses were again sourced from Audi & Wapstra (1995). In practice, the ratio of the square root terms enters the opacity anyway through the mass scalefactor, and the microturbulent velocity is not used in 3D (i.e.  $\xi_{\rm t}=0$ was used) so these can be disregarded also. The total opacity scalefactor applied in each case was therefore the product of the appropriate one of the reduced mass correction factors and the ratio of reference (12C or 16O) isotope to that involved in the line being modelled (13C or 18O), i.e.

\varsigma_{\kappa} = \frac{N(^b{\rm X})}{N(^a{\rm X})}\left...
\end{displaymath} (B.12)

Had the opacity scalefactor reflected abundance differences only, not intrinsic atomic (i.e. reduced mass) corrections, our adopted 3D-based ratios would have been 12C/13C = 92.9 +4.1-3.9 and 16O/18O = 516 +32-30.

The opacity scalefactors were then iteratively altered in calculations of the 13C16O and 12C18O line lists to reflect different isotopic ratios, in the same manner as the bulk carbon abundance was iterated to determine the "12C'' abundance using the three sets of 12C16O lines. The same bulk carbon abundance as used in the final iteration of the weak 12C16C calculation was used in each case, and the isotopic ratios after some iteration i given by

\left(\frac{^a{\rm X}}{^b{\rm X}}\right)_{i} = \left(\frac{...
...lta ^{12}{\rm C}_{\rm final}}}{10^{\Delta {\rm iso}_{i}}}\cdot
\end{displaymath} (B.13)

Here $\Delta {\rm iso}_{i}$ and $\Delta ^{12}{\rm C}_{\rm final}$ are the abundance corrections produced by the ith isotopomer iteration and the final iteration of whichever 12C16O list is used to indicate the 12C abundance, respectively.

B.3 Fractional scalefactor

Following opacity scalefactor convergence, a further iteration was also performed on each of the sets of 12C16O lines using a fractional scalefactor. This accounts for the fact that not all carbon in CO lines exists in 12C16O. This scalefactor was calculated from the just derived 12C/13C and 16O/18O ratios, such that

\begin{eqnarray*}\varsigma_{\rm frac} = \frac{N(^{12}{\rm C}^{16}{\rm O})}{N({\rm CO})},

and since

\begin{eqnarray*}N({\rm CO}) = N(^{12}{\rm C}^{16}{\rm O})+N(^{13}{\rm C}^{16}{\rm O})+N(^{12}{\rm C}^{18}{\rm O})+N(^{12}{\rm C}^{17}{\rm O})

we see that
                           $\displaystyle %
\varsigma_{\rm frac}$ = $\displaystyle \left(1 + \frac{N(^{13}{\rm C}^{16}{\rm O})}{N(^{12}{\rm C}^{16}{\rm O})} + \frac{N(^{12}{\rm C}^{18}{\rm O})}{N(^{12}{\rm C}^{16}{\rm O})} \right.$  
    $\displaystyle \left. + \frac{N(^{12}{\rm C}^{17}{\rm O})}{N(^{12}{\rm C}^{16}{\rm O})}\right)^{-1}$  
  = $\displaystyle (1 + {^{13}{\rm C}}/{^{12}{\rm C}} + {^{18}{\rm O}}/{^{16}{\rm O}} + {^{17}{\rm O}}/{^{16}{\rm O}})^{-1}.$ (B.14)

Therefore, $\varsigma_{\rm frac}$ represents the fraction of the bulk carbon abundance actually indicated by 12C16O lines. Seeing as the 12C17O lines are so weak in the ATMOS spectrum, accurate derivation of the 16O/17O ratio would not have been possible in this study, so the value of 16O/17O used in the calculation of fractional scalefactors was the terrestrial value of $\sim$2630 (Coplen et al. 2002). Being so incredibly small, the contribution of any difference in the 17O abundance between Earth and the Sun would have had almost no effect upon the resultant scalefactor. Using these fractional scalefactors, the final iterations of the 12C16O line lists indicated bulk solar carbon abundances.

Appendix C: Line lists

Refer to Tables C.1-C.6.

Table C.1: Strong 12C16O line list used for the Sect. 5 study, consisting of 31 lines. $\log gf$ data calculated from Goorvitch & Chackerian (1994).

Table C.2: Weak 12C16O line list used for the Sect. 6 study, consisting of 13 lines. $\log gf$ data calculated from Goorvitch & Chackerian (1994).

Table C.3: 13C16O line list used for the Sect. 6 study, consisting of 16 lines. $\log gf$ data calculated from Goorvitch & Chackerian (1994).

Table C.4: 12C18O line list used for the Sect. 6 study, consisting of 15 lines. $\log gf$ data calculated from Goorvitch & Chackerian (1994).

Table C.5: Low excitation (LE) 12C16O line list used for the Sect. 6 study, consisting of 15 lines. $\log gf$ data calculated from Goorvitch & Chackerian (1994).

Table C.6: First overtone ( $\Delta v = 2$) 12C16O line list used for the Sect. 6 study, consisting of 66 lines. $\log gf$ data calculated from Goorvitch & Chackerian (1994).

Copyright ESO 2006