Two-dimensional non-LTE \ion{O}{I} 777\,nm line formation in radiation hydrodynamics simulations of Cepheid atmospheres

TL;DR

This study assesses the accuracy of 1D hydrostatic models for oxygen abundance measurements in Cepheids by comparing them with more realistic 2D non-LTE radiation hydrodynamics simulations, revealing significant phase-dependent differences.

Contribution

It demonstrates the limitations of 1D LTE models for oxygen line analysis in Cepheids and highlights the importance of phase selection for accurate abundance determination.

Findings

2D non-LTE vs 1D LTE abundance differences range from -1.0 to -0.25 dex.

2D non-LTE vs 1D non-LTE differences range from -0.2 to 0.8 dex.

Best phases for oxygen abundance analysis are between 0.3 and 0.8, near hydrostatic equilibrium.

Abstract

Oxygen abundance measurements are important for understanding stellar structure and evolution. Measured in Cepheids, they further provide clues on the metallicity gradient and chemo-dynamical evolution in the Galaxy. However, most of the abundance analyses of Cepheids to date have been based on one-dimensional (1D) hydrostatic model atmospheres. Here, we test the validity of this approach for the key oxygen abundance diagnostic, the \ion{O}{I} $777\,\mathrm{nm}$~triplet lines. We carry out 2D non-LTE radiative transfer clculations across two different 2D radiation hydrodynamics simulations of Cepheid atmospheres, having stellar parameters of $T_\mathrm{eff}= 5600$ K, solar chemical compositions, and $\log\,g= 1.5$ and $2.0$, corresponding to pulsation periods of 9 and 3 days, respectively. We find that the 2D non-LTE versus 1D LTE abundance differences range from $-1.0$~dex to $-0.25$~dex depending on pulsational phase. The 2D non-LTE versus 1D non-LTE abundance differences range from $-0.2$~dex to $0.8$~dex. The abundance differences are smallest when the Cepheid atmospheres are closest to hydrostatic equilibrium, corresponding to phases of around $0.3$ to $0.8$, and we recommend these phases for observers deriving the oxygen abundance from \ion{O}{I} $777\,\mathrm{nm}$ triplet with 1D hydrostatic models.

Full text

1 Introduction

Oxygen is the third most abundant chemical element and the most abundant metal in the Universe. As with other $\upalpha$-elements, it forms primarily in massive stars. As such, oxygen abundances in the atmospheres of stars are important tracers of the Galactic chemical evolution (GCE; e.g. Chiappini et al., 2003; Minchev et al., 2013; Kubryk et al., 2015a, b). In particular, they provide strong constraints on the formation timescales of Galactic stellar populations (for details, see e.g. Ramírez et al. 2007, 2013 and Magrini et al. 2017).

In this context, classical Cepheids are interesting; as they are young (¡250 Myr; e.g. Bono et al., 2005), massive, evolved stars. Takeda & Takada-Hidai (1998) and later Takeda et al. (2013) investigated the nature of evolution-induced mixing in the envelope of evolved intermediate-mass stars in supergiants and Cepheid variables of various pulsation periods to determine the photospheric abundances of C, N, O, and Na, applying non-LTE analyses. They showed that the observed CNO abundance trends are mainly a result of canonical dredge-up of CN-cycled material, while any significant non-canonical deep mixing of ON-cycled gas is ruled out.

With the exception of a few elements whose composition is modified by the first dredge-up (Li, C, and N in particular), the chemical composition of Cepheids reflects the current chemical composition of the interstellar medium (ISM). Furthermore Cepheids are bright stars that can be detected out to large distances. They can be used as standard candles for distance measurements due to the existence of a relation between their period of pulsation and intrinsic luminosity (Leavitt, 1908; Leavitt & Pickering, 1912). Combining information on the oxygen abundance in Cepheids with their distances, we can derive the present-day oxygen abundance gradient of the Galactic disc, which is a crucial constraint for the models of the Galactic chemical evolution. In this respect, Cepheids adequately complement large-scale surveys such as GALAH (De Silva et al., 2015), APOGEE (Majewski et al., 2017), and GAIA-ESO (Gilmore et al., 2012), focusing mostly on older stars; they are therefore useful for deriving the time evolution of Milky Way gradients (e.g. Anders et al., 2017).

A number of studies were performed (e.g. Lemasle et al., 2008, 2013; Luck et al., 2011; Luck & Lambert, 2011; Genovali et al., 2015; da Silva et al., 2016) to determine Galactic abundance gradients for numerous $\upalpha$-, iron-peak, and neutron-capture elements using a large number of classical Cepheids. However, the abundance measurements presented in these Cepheid studies are all based on 1D hydrostatic stellar model atmospheres. In these models, convection is typically treated in the framework of the mixing length theory (Böhm-Vitense, 1958). Abundances are typically measured at different phases, and thereby represent the dynamical, evolving atmosphere as a sequence of quasi-hydrostatic states with different effective temperatures and surface gravities. Furthermore, in most of the above-mentioned studies, the line formation was assumed to take place under the condition of local thermodynamic equilibrium (LTE), although later studies relaxed this assumption and were based on non-LTE radiative transfer.

Possible non-LTE effects in Cepheids were discussed by Luck & Lambert (2011). Later, Luck et al. (2013) and Korotin et al. (2014) presented 1D non-LTE oxygen abundance analyses of large samples of Cepheids. The O I $777\,\mathrm{nm}$ and $845\,\mathrm{nm}$ triplets can be used to derive oxygen abundances because they are easily detectable over a wide range of effective temperature and generally do not suffer from severe blending with other species; furthermore, the spectral region is less contaminated by terrestrial bands compared to the forbidden low-excitation [O I] $630.0\,\mathrm{nm}$ and $636.4\,\mathrm{nm}$ lines, or to the permitted high-excitation O I $615.8\,\mathrm{nm}$ triplet. However, the near-IR triplets are sensitive to departures from LTE. For a given oxygen abundance, taking into account non-LTE effects tends to make these lines stronger; thus, the inferred oxygen abundance is lower once the assumption of LTE is relaxed.

Since calculations of multi-dimensional Cepheid models are a sizable computational problem due to the presence of the different spatial and temporal scales, it is not a surprise that the application of 1D hydrostatic stellar model atmospheres is common. However, it is desirable to obtain an understanding of the limitations inherent to the standard approach. Working towards this goal, in a previous paper (Vasilyev et al., 2017) we presented a two-dimensional (2D) model of a short-periodic Cepheid-like variable. The model has realistic convection based on the conservation laws. Additionally, the model shows self-excited fundamental mode pulsations due to the $\kappa$-mechanism (Eddington, 1917; Zhevakin, 1963). In the present paper, performing similar 2D radiation-hydrodynamics (RHD) simulations we aim to investigate non-LTE effects on oxygen, their phase-dependence, and biases in the derived oxygen abundances compared to a standard 1D non-LTE analysis based on hydrostatic models.

Specifically, in Section 2 we describe the line formation code, atomic model, as well as 1D and 2D model atmospheres. This study focuses on one of the key oxygen abundance diagnostics, namely the O I $777\,\mathrm{nm}$ triplet. The results including line profiles, abundance errors and biases of the standard approach are presented in Section 3. In Section 4 we discuss the optimal phase for observers in order to get unbiased measurements of the oxygen abundance in Cepheids. We also compare the behaviour of departures from LTE in our 2D and 1D model atmospheres.

2 Methods

2.1 2D RHD simulations of Cepheid variables

We performed 2D RHD simulations of Cepheid-like variables using the RHD code CO5BOLD (Freytag et al., 2012). The code solves the set of 2D compressible RHD equations using a short-characteristics method for the radiative transfer, which was considered in grey approximation. Grey PHOENIX (Allard et al., 2012) merged with OPAL (Iglesias et al., 1992) opacities were used in the simulations. In this study two models are considered: a three-day model (M3) and a nine-day model (M9). The models have a mean effective temperature $T_{\mathrm{eff}}=5600\,\mathrm{K}$, constant depth-independent gravitational acceleration, $\log\,g=1.5$ for M9 and $\log\,g=2.0$ for M3, and solar-like chemical composition. In Vasilyev et al. (2017) we discussed the motivation and limitations of our approach.

The model M3 has a three-day period of pulsation; it was introduced and analysed in Vasilyev et al. (2017). In the present study we extended the sequence to $14$ full pulsation periods, compared to 6 in our previous paper. We illustrate the bolometric light curve of the model in the left panel of Fig. 1. Due to the impact of the stochastic convective component, the light curve shows cycle-to-cycle variations (for details, see Vasilyev et al., 2017).

After relaxation of initial perturbations, model M9 reaches an approximate limit cycle exhibiting fundamental-mode oscillations with a period of approximately $9$ days. We illustrate the bolometric light curve of the model in the right panel of Fig. 1. Due to the impact of convective noise, the effective temperature of the phases of maximum light varies between $6460\,\mathrm{K}$ and $6730\,\mathrm{K}$; in particular, during the seventh cycle the maximum effective temperature reaches a value above $6700\,\mathrm{K}$.

While the 2D models have a depth-independent gravitational acceleration $\@vec{g_{\mathrm{o}}}$, spectral lines are in fact sensitive to the effective gravitational acceleration, $\@vec{g}_{\mathrm{eff}}=\@vec{g_{\mathrm{o}}}-\frac{d\@vec{v}}{dt}$, which is balanced by the pressure gradient. For the maximum of the seventh cycle of the M9 model, the effective gravitational acceleration in the photosphere drops by a factor of around three due to the impact of the convective noise. This produces a loop-like structure in our results for the abundance corrections and abundance errors, as we discuss in detail in Sect. 3.

Both simulations use an equidistant spatial Cartesian grid of $N_{x}\times N_{z}=600\times 500$ along the horizontal and vertical directions. The geometrical extensions of the models along the horizontal and vertical directions are $l_{\mathrm{x}}\times l_{\mathrm{z}}=1.1\cdot 10^{12}\times 0.5\cdot 10^{12}$ cm and $l_{\mathrm{x}}\times l_{\mathrm{z}}=3.6\cdot 10^{12}\times 1.7\cdot 10^{12}$ cm for M3 and M9, respectively. The models both have periodic horizontal boundary conditions, and closed top and bottom boundaries. The pulsations induce strong perturbations in the vicinity of the upper (closed) boundary. In order to damp these perturbations a drag force was applied on the horizontal mean motion in the uppermost few grid layers. The impact of this approach on spectroscopic properties is discussed in Vasilyev et al. (2017).

The effective temperature of the 2D models $T_{\mathrm{eff}}^{\mathrm{2D}}$ was calculated from the bolometric radiative flux passing the top boundary at each instance in time. The kinematic acceleration $\frac{d\@vec{v}}{dt}$, which is variable during the pulsational cycle, leads to a variation in $\@vec{g}_{\mathrm{eff}}$ in the range $\approx(0.6\ldots 3.5)\@vec{g_{\mathrm{o}}}$ for M9, and $\approx(0.5\ldots 3.0)\@vec{g_{\mathrm{o}}}$  for M3. It was calculated from the motion of a (Lagrangian) mass element located close to optical depth unity. We illustrate the variation in $T_{\mathrm{eff}}$ and $\log\,g_{\mathrm{eff}}$ over the pulsational cycles for models M3 and M9 in Fig. 2.

2.2 One-dimensional models

We constructed a grid of 1D hydrostatic model stellar atmospheres using the Lagrangian hydrodynamics (LHD) code, which solves the set of 1D RHD equations in the Lagrangian frame. In the present project we operated it in hydrostatic mode, which provides structures in hydrostatic equilibrium. To facilitate a differential comparison with the 2D models, the same opacities and equation-of-state tables were used as in the 2D models. The dimensionless mixing-length parameter $\alpha$ was fixed to $1.5$ for all 1D models. In order to cover the parameters encountered in M3 and M9, $T_{\mathrm{eff}}$, and $\log\,g$ for the 1D models were chosen in the range $5200<T_{\mathrm{eff}}/\mathrm{K}<6800$ with a step $\Delta T_{\mathrm{eff}}=100\,\mathrm{K}$ and $1.1<\log\,g<2.9$ with $\Delta\log\,g=0.2$, respectively. We depict the grid of 1D models in the $\log\,g$—$T_{\mathrm{eff}}$ plane in Fig. 2 which also shows the parameters for the model M3 and M9. Technical details concerning 1D hydrostatic models can also be found in Vasilyev et al. (2018).

2.3 Microturbulent velocity

To calculate 2D LTE/non-LTE versus 1D LTE/non-LTE abundance corrections, it is necessary to first consider the microturbulent velocity that needs to be adopted in the spectral synthesis based on 1D model atmospheres. The microturbulent velocity is a relevant parameter in spectroscopic analyses of late-type stars (both pulsating and non-pulsating) that are based on 1D hydrostatic models. It is needed to account for the additional line broadening induced by gradients in the (convective) velocity fields over distances $\Delta\tau_{\mathrm{R}}\lesssim 1$ (for details, see e.g. Steffen et al., 2013). Due to the change in the thermal structure and velocity gradients during the pulsational cycle, the microturbulent velocity is phase-dependent (for details, see Luck & Andrievsky 2004, Andrievsky et al. 2005, and Kovtyukh et al. 2005, and the later studies of Vasilyev et al. 2017 and Proxauf et al. 2018a).

For each 2D snapshot, we calculated an effective microturbulent velocity, which we used in the 1D line formation calculations (Sect. 2.4) for a comparison of 2D and 1D results. Following Steffen et al. (2013) and Vasilyev et al. (2017), these microturbulent velocities were inferred from 2D (LTE/non-LTE) spectra. The microturbulent velocity in the LTE case is calculated by equating equivalent widths $W(\@vec{v}_{\mathrm{2D}})=W(\xi_{\mathrm{t}})$, where $W(\@vec{v}_{\mathrm{2D}})$ is calculated with the full 2D velocity field $\@vec{v}_{\mathrm{2D}}$, and $W(\xi_{\mathrm{t}})$ by replacing $\@vec{v}_{\mathrm{2D}}$ by a depth-independent isotropic microturbulent velocity $\xi_{\mathrm{t}}$.

As in the case of observations, Fe i and Fe ii lines were used for the diagnostic of the microturbulent velocity. The line list, which was given in Vasilyev et al. (2018), contains artificial neutral and singly ionised iron lines with excitation potentials $1$, $3$, and $5\,\mathrm{eV}$. Oscillator strengths were chosen to cover the linear and non-linear parts of the curve of growth. The spectral synthesis was performed using the code LINFOR3D (Gallagher et al., 2017). Measured microturbulent velocities were averaged over line strengths and excitation potentials within ionisation specie for each individual 2D snapshot. Details concerning the phase-dependence of microturbulent velocities of Fe i and Fe ii lines as a function of the photometric phase in Cepheids will be presented in forthcoming papers.

2.4 Atomic model and non-LTE line formation

Amarsi et al. (2018a) presented a new oxygen model atom that incorporates inelastic O+H collisional rate coefficients based on the asymptotic two-electron model of Barklem (2018) and the free electron model of Kaulakys (1991). The MPI-parallelised 3D non-LTE radiative transfer code Balder (Amarsi et al., 2018b), a custom version of MULTI3D (Leenaarts & Carlsson, 2009), was used in this study. The code solves the statistical equilibrium equation taking into account radiative and collisional bound-bound and bound-free transitions to get level populations and also to calculate the final emergent spectrum. Line formation calculations were performed on $705$ snapshots of M9, spanning $10$ cycles, and $326$ snapshots of M3, spanning 14 cycles.

The LTE and non-LTE calculations on the snapshots of the 2D models (Sect. 2.1) and the LTE and non-LTE calculations across the grid of 1D models (Sect. 2.2) proceeded in mainly the same way. For the 1D calculations, a depth-independent microturbulent velocity $\xi_{\mathrm{t}}$ was introduced, and was varied between $0.0$ and $4.0\,\mathrm{km\,s^{-1}}$ with a step $1.0\,\mathrm{km\,s^{-1}}$. Calculations were performed for several oxygen abundances $-0.8<\mathrm{[O/Fe]}<+0.8$ with a step $0.4\,\mathrm{dex}$.

2.5 Abundance corrections

In order to compare the analysis using multi-dimensional model atmospheres with the classical 1D approach we can calculate abundance corrections and abundance errors. Given a measured LTE abundance, we define the abundance correction as the abundance difference compared to LTE, {IEEEeqnarray}rCl Δ^NLTE_LTE (log  ϵ^LTE)&= log  ϵ^NLTE - log  ϵ^LTE ,

such that the non-LTE equivalent width matches the LTE equivalent width.

3 Results

3.1 Departure coefficients

In order to understand the behaviour of the abundance corrections $\Delta^{\mathrm{2D\,NLTE}}_{\mathrm{1D\,NLTE}}$ and $\Delta^{\mathrm{2D\,NLTE}}_{\mathrm{1D\,LTE}}$, we can look at the variation of the departure coefficients, $b(\tau_{\mathrm{R}})=n_{\mathrm{NLTE}}/n_{\mathrm{LTE}}$, for the 1D hydrostatic models and 2D models. For the first, $b^{\mathrm{1D}}$ values were calculated for the lower $3\mathrm{s}^{5}\,\mathrm{S}$ and upper $3\mathrm{p}^{5}\,\mathrm{P}$ levels of the O i 777 nm triplet in the 4D parameter space $\tau_{\mathrm{R}},T_{\mathrm{eff}},\log\,g_{\mathrm{eff}},\xi_{\mathrm{t}}$. For each combination \Big{[}T_{\mathrm{eff}},\log\,g_{\mathrm{eff}},\xi_{\mathrm{t}}\Big{]}^{\mathrm{2D}} a cubic interpolation that splits the multi-dimensional case into a sequence of 1D interpolations was performed on the grid of 1D models to obtain departures from LTE as functions of the Rosseland optical depth.

In the 2D models – due to the horizontal inhomogeneities produced by convection – departures from LTE $b(\tau_{\mathrm{R}})$ are different for in- and outflows (e.g. Amarsi et al., 2016). In order to compare departures from LTE for 2D models with the 1D case, $b^{\mathrm{2D}}$ values for each individual 2D snapshot were horizontally averaged over surfaces of constant optical depth. Both 2D models predict a similar dependence of the abundance corrections on effective temperature. Since they are similar we focused on the nine-day model, and calculated median 2D departures from LTE.

We plot the departure coefficients in Fig. 3. Namely, the median 2D departures from LTE as a function of $\tau_{\mathrm{R}}$ and photometric phase for the lower and upper levels are shown by solid blue and red lines, respectively. One-dimensional interpolated $b^{\mathrm{1D}}$ departures from LTE are shown by dashed lines using the same colours for the levels. Qualitatively, the departures from LTE of 1D and 2D models show a similar dependence on the photometric phase. The lower and upper levels of the O I $777\,\mathrm{nm}$ triplet, namely $3\mathrm{s}^{5}\,\mathrm{S}$ and $3\mathrm{p}^{5}\,\mathrm{P}$, are almost always overpopulated during the pulsating cycle. At photometric phases $\phi_{\mathrm{ph}}\approx 0.3\ldots 0.8$ the lower and upper levels are less underpopulated compared to the rest of the phases, due to the minimum on the light curve and effective temperature. Quantitatively, for this range of phases 1D and 2D models predict very similar optical depth dependences of the departures from LTE. As a result, the abundance correction $\Delta^{\mathrm{2D\,NLTE}}_{\mathrm{1D\,NLTE}}$ is small. Moreover, the 2D non-LTE versus 1D LTE abundance correction turns out to be smallest, as shown in Sect. 3.5.

3.2 O I $777\,\mathrm{nm}$ triplet profiles

We illustrate the temporal evolution (over a pulsation cycle) of the normalised emergent LTE and non-LTE O I $777\,\mathrm{nm}$ triplet flux profiles in Fig. 4 for M3 and M9. The line profiles are shown without adding instrumental broadening or noise. This allows us to identify the occasional multi-component structure of the line profiles, for example at photometric phases $\phi_{\mathrm{ph}}=0.299$ and $0.318$ for M3 and M9, respectively. Overall, expansion or compression produce asymmetric line profiles that are typical for pulsating stars (Shapley & Nicholson, 1919).

In Fig. 5 we show the equivalent widths as a function of the effective temperature and the surface gravity, for both M3 and M9 models. Generally, the non-LTE equivalent widths are larger than the LTE values, as also seen in Fig. 4. The main driver of the non-LTE effect is photon losses in the O I $777$ nm triplet itself (e.g. Takeda, 2003), which leads to a strengthening of the triplet.

Figure 5 also shows that the equivalent widths have a strong dependence on the effective temperatures, and a weaker dependence on the surface gravity. The models span a wide range of effective temperatures, and a relatively narrow range of surface gravities (Fig. 2): at fixed $\log\,g$, $T_{\mathrm{eff}}$ spans up to around $1000\,\mathrm{K}$, corresponding to around a $0.3$ to $0.5\,\mathrm{dex}$ change in line strength; instead, at fixed $T_{\mathrm{eff}}$, $\log\,g$ spans up to around $0.4\,\mathrm{dex}$, corresponding to a less than $0.04\,\mathrm{dex}$ change in line strength, based on the 1D model atmospheres. In other words, in this context the line strengths are more sensitive to the effective temperature than to the surface gravity, as expected for lines of high excitation potential from a neutral species.

3.3 Abundance differences: 2D non-LTE versus 2D LTE

In Fig. 6 we illustrate the 2D LTE versus 2D non-LTE abundance corrections as functions of effective temperature and surface gravity for the 2D models. There is a larger spread in the abundance corrections at fixed surface gravity than at fixed effective temperature. Furthermore the M3 and M9 models are clearly offset in abundance corrections versus surface gravity space, whereas they overlap in in abundance corrections versus effective temperature space. This can be understood from our discussion in Sect. 3.2: the non-LTE effect is one of photon losses in the lines themselves, and the triplet is more sensitive to the effective temperature than to the surface gravity.

Figure 6 shows that the non-LTE effects are the most severe at higher effective temperatures. This happens close to the phases of maximum light; as seen in Fig. 4, at these phases the non-LTE line depths are roughly twice as large as the LTE values. The M3 model has smaller line depths ratios between the LTE and non-LTE cases because at maximum light it has a smaller effective temperature (and thus weaker O I $777$ nm lines and fewer photon losses). At minimum light (low effective temperature), the M3 and M9 models have similar line depth ratios in both LTE and non-LTE, with less severe non-LTE effects. Hence, there is a modulation of the non-LTE effects over the pulsations. Nevertheless, even at the phase of maximum expansion and start of compression (minimum of the effective temperature) the 2D LTE versus 2D non-LTE abundance errors are of the order of $+0.6$ dex, and become larger up to $+1.3$ dex at $T_{\mathrm{eff}}\approx 6700$ K. In order to obtain unbiased estimates of the oxygen abundance in Cepheids (based on the O I $777$ nm triplet), we have to consider the effects of departures from LTE.

The scatter of the abundance correction at fixed effective temperature is related to the impact of the convective noise. It is expected to become weaker when an average over a larger number of convective elements is performed, for example in 3D models or by averaging across many cycles. The M3 model shows a larger scatter due to the stronger impact of the convective component at higher surface gravity.

3.4 Abundance differences: 1D non-LTE versus 1D LTE

We illustrate variations in the 1D non-LTE versus 1D LTE abundance corrections for the $777.4\,\mathrm{nm}$ component of the O I $777\,\mathrm{nm}$ triplet as a function of effective temperature for 1D models in Fig. 7. Luck et al. (2013) derived abundance corrections for an ensemble of Galactic Cepheids for the O i $777.4\,\mathrm{nm}$ line as function of effective temperature which are shown as black stars in Fig. 7. As can be seen, their results are qualitatively similar to ours; however, the temperature dependence of their abundance corrections is steeper.

This could be due to several factors. First of all, M3 and M9 are not complete stellar models. Their pulsational amplitudes and periods depend to some extent on the applied geometry and boundary conditions. This leads to a different $T_{\mathrm{eff}}-\log\,g-\mathrm{Amplitude}-\mathrm{Period}$ relation with respect to what is expected for real Cepheids. Nevertheless, we note here that when focusing on spectroscopic properties certain observed spectroscopic features can be reproduced by the 2D models (Vasilyev et al., 2017).

Secondly, the surface gravity of our 1D models is based on the effective acceleration $\@vec{g}_{\mathrm{eff}}$ of an Lagrangian mass element located near optical depth unity in the 2D models. This differs from a $\@vec{g}_{\mathrm{eff}}$ adopting the condition of ionisation balance, which is traditionally used in classical 1D analyses to measure the surface gravity (Vasilyev et al., 2018).

Thirdly, Luck et al. (2013) used 1D hydrostatic ATLAS9 model atmospheres (Kurucz, 1992), whereas we used the 1D LHD model atmospheres described in Sect. 2.2; the use of grey radiative transfer has a significant impact on the temperature stratification of the LHD models, and it is possible that this may be the reason for most of the differences. Differences in the atomic models may also contribute.

Additionally, the result of Luck et al. (2013) shows a larger vertical scatter for fixed effective temperature. They explain it by the sensitivity of the oxygen abundance to the microturbulent velocity and strengths of the line.

3.5 Abundance differences: 2D versus 1D

In this section we compare the 2D and 1D approaches to determine the optimal phases for measuring the O I $777\,\mathrm{nm}$ triplet such that the 2D versus 1D differences are at a minimum. To associate a 1D model with a given 2D snapshot, we adopt an effective temperature and effective surface gravity as described in Sect. 2.1, and a microturbulent velocity as described in Sect. 2.3. As stated before, this does not exactly correspond to the standard analysis procedure where the effective temperature is typically determined from line depth ratios (Gray, 1994; Kovtyukh & Gorlova, 2000), and the gravity from the ionisation balance of Fe i to Fe ii; the microturbulent velocity is obtained from the condition of a consistent (often iron) abundance derived for weak and strong lines (in LTE).

In our case, for each 2D snapshot the equivalent widths from the grid of 1D models were interpolated for the 2D parameters $\left[T_{\mathrm{eff}},\log\,g_{\mathrm{eff}},\xi_{\mathrm{t}}\right]^{\mathrm{2D}}$, and the abundance corrections and abundance errors were then calculated as stated earlier. The multi-dimensional interpolation was split into a sequence of 1D cubic interpolations. For the 2D non-LTE versus 1D non-LTE abundance corrections we illustrate the results of the interpolation for M3 and M9 in the upper panel of Fig. 8. The total equivalent widths of the O I $777\,\mathrm{nm}$ triplet were taken in the calculations here. Analogously, 2D non-LTE versus 1D LTE abundance corrections are presented in the lower panel of Fig. 8.

In Fig. 8 models M3 and M9 show similar slopes for the 2D non-LTE versus 1D non-LTE abundance corrections as functions of effective temperature. The slope does not depend on the surface gravity and is primarily controlled by the change of the thermal structure over the pulsational cycle. Interestingly, we find that at certain phases a 1D non-LTE analysis would provide the same abundance as a 2D non-LTE analysis. Qualitatively speaking, these phases are near maximum expansion and start of compression. From these phases, at higher effective temperatures, the bias grows such that the 1D non-LTE analysis gives abundances that are as much as $0.5\,\mathrm{dex}$ too low.

In Fig. 8 models M3 and M9 show similar slopes for the 2D non-LTE versus 1D LTE abundance corrections as functions of effective temperature. As with the 2D non-LTE versus 1D non-LTE abundance corrections, the slope does not depend on the surface gravity and is primarily controlled by the change in the thermal structure over the pulsational cycle. Unlike the situation for the 2D NLTE versus 1D NLTE corrections, however, we do not find any phases at which a 1D LTE analysis would provide the same abundance as a 2D non-LTE analysis. In the best case, at certain photometric phases, corresponding to effective temperatures $T_{\mathrm{eff}}\lesssim 5700\,\mathrm{K}$, the abundance corrections are almost insensitive to the change in $T_{\mathrm{eff}}$ and are around $-0.3\,\mathrm{dex}$, ranging from $-0.2\,\mathrm{dex}$ to $-0.4\,\mathrm{dex}$. For higher effective temperatures the abundance corrections become more severe, decreasing to $-1.5\,\mathrm{dex}$.

In Fig. 8, for certain phases the abundance corrections show loop-like structures. As we discussed in Sect. 2.1, for one cycle at maximum light, $\@vec{g}_{\mathrm{eff}}$ drops by a factor of three while ${T}_{\mathrm{eff}}$ increases by $130$ K. For this short time interval, the increase in the O I $777\,\mathrm{nm}$ triplet component equivalent widths in LTE and non-LTE cases are primarily due to the sensitivity to the changes in $\log\,g_{\mathrm{eff}}$. In our calculations, the relative changes of LTE and non-LTE equivalent widths are $\left(\frac{\Delta W}{W}\right)_{\mathrm{LTE}}\approx 0.4$ and $\left(\frac{\Delta W}{W}\right)_{\mathrm{NLTE}}\approx 0.2$, respectively. As a result, the 2D non-LTE versus 1D LTE abundance error increases from $-1.35\,\mathrm{dex}$ to $-1.1\,\mathrm{dex}$.

4 Discussion

4.1 Optimal phases for 1D analyses

In practice, spectroscopic analyses of the O I $777\,\mathrm{nm}$ triplet are based on 1D hydrostatic model atmospheres, and will continue to be for the foreseeable future. However, based on the 2D non-LTE abundance corrections compared to 1D non-LTE and 1D LTE, we can provide recommendations on which phases to target, to obtain the most reliable oxygen abundances. Our recommendation for 1D LTE and 1D non-LTE analysis is to avoid phases around maximum light as shown in Fig. 9. In a 1D LTE analysis, in the best case we can obtain abundance corrections of the order of $-0.3\pm 0.1\,\mathrm{dex},$ and in the case of 1D non-LTE analysis $0\pm 0.1\,\mathrm{dex}$ at photometric phases $\phi_{\mathrm{ph}}\approx 0.3\ldots 0.8$ independent of period. Within these phases, disturbances of the thermal structure and convection are not as strong as during maximum compression. As a result, the atmosphere is roughly in hydrostatic equilibrium, and the 1D models reasonably reproduce the mean thermal structure of the line formation region of the 2D model. If possible, a 1D non-LTE analysis based on these phases would be preferable.

4.2 Implications on Galactic oxygen abundances

In order to investigate the impact of the phase variation of the 2D non-LTE versus 1D non-LTE abundance corrections on the Galactic oxygen abundance gradient, 1D non-LTE oxygen abundances for Galactic Cepheids based on the O I 777 nm triplet were taken from Korotin et al. (2014). Using a linear fit we measured the gradient for two cases: (i) spectra were taken at all phases and (ii) at optimal phases $0.3<\phi_{\mathrm{ph}}<0.8$ for which a 1D non-LTE analysis should provide largely unbiased abundances according to our predictions. The gradients are presented in Fig. 10 as blue and red lines, respectively, including $1\sigma$ uncertainties. For the data obtained at all phases the slope is $-0.053\pm 0.003$ dex$\cdot$kpc*-1*, whereas with data at optimal phases we get a slightly steeper gradient $-0.056\pm 0.004$ dex$\cdot$kpc*-1*. However, the change in slope is within the uncertainties.

It is also possible to apply the 2D non-LTE versus 1D non-LTE abundance corrections directly to a 1D analysis; however, this requires 2D models with different oxygen abundances and metallicites, as well as pulsation properties, in order to cover the observed range of parameters. Recent multi-phase analyses of Galactic Cepheids by Luck (2018) and Proxauf et al. (2018b) using 1D model atmospheres show some dependence of elemental abundances on the photometric phase. In particular, a dependence is seen in the iron abundances, and they are possibly an indication of the failure of 1D hydrostatic models to represent the 3D hydrodynamic nature of a pulsating Cepheid atmosphere. For oxygen, the phase dependence should be stronger for the the O I 777 nm triplet owing to its temperature sensitivity, which is amplified by temperature-sensitive 3D non-LTE effects. We expect the dependences to be stronger in Cepheids with higher mean effective temperatures, to be more metal-rich, and to have longer periods.

5 Conclusion

In this paper for the first time we checked the validity of the standard non-LTE spectroscopic analysis of oxygen in the context of Cepheid variables. To achieve this goal, we carried out 2D RHD simulations pulsating Cepheid atmospheres with three- and nine-day periods. In post-processing, we performed 2D non-LTE radiative transfer calculations for atomic oxygen. We compared the 2D results to those based on 1D hydrostatic model atmospheres. We find a strong modulation of the non-LTE effects in the oxygen 777 nm triplet lines during the pulsations phase, mostly due to the change of the effective temperature. The dependence of the abundance differences $\Delta^{\mathrm{2D\,NLTE}}_{\mathrm{2D\,LTE}}$ based on the 2D Cepheid models is different with respect to the analysis of Luck et al. (2013) obtained with the non-LTE 1D hydrostatic approach. We confirmed that 1D LTE analyses of the O I 777 nm triplet cannot provide unbiased estimates of the oxygen abundance at any phase. This also implies that 2D effects do not cancel out the non-LTE effects.

The 1D non-LTE analysis of the O i 777 nm triplet lines can provide an unbiased estimate of the oxygen abundance only for a limited range of photometric phases. For the analysis using the whole oxygen 777 nm triplet, it is desirable to analyse spectra taken at photometric phases $\phi_{\mathrm{ph}}\approx 0.3\ldots 0.8,$ independent of period, when the pulsating atmosphere is closer to hydrostatic conditions. For these phases 1D as well as 2D non-LTE calculations predict a similar dependence of the departure coefficients for the lower and upper level on the optical depth.

In this work we presented a first step towards more reliable spectroscopic analyses of Cepheid variables. Nevertheless, there are a number of ways in which the models still have to be improved. The model atmospheres themselves are 2D: for non-variable stars there is evidence that full 3D simulations are required to generate the correct velocity fields and temperature and pressure stratifications (Asplund et al., 2000); we expect that similar problems may exist here. In particular, with full 3D simulations we would expect less influence of the stochastic convective noise on pulsating properties and thermal structure of the atmosphere. It reduces cycle-to-cycle variations of the light curve and ultimately scatter, which we see near the maximum light phases in LTE versus non-LTE abundance corrections at fixed effective temperature. Additional avenues of work involve constructing the atmospheres using non-grey radiative transfer, perhaps by employing opacity bins. Another important step is to broaden the model basis to better sample the $T_{\mathrm{eff}}$ – $\log\,g$ – [Fe/H] – Amplitude – Period space. This is a necessary step, before 2D non-LTE (or 3D non-LTE) can be directly applied to the literature 1D LTE abundances.

Acknowledgements.

VV, HGL, and BL acknowledge financial support from the Sonderforschungsbereich SFB 881 ‘The Milky Way System’ (subprojects A4, A5) of the German Research Foundation (DFG) and MPIA IT department for computer resources. AMA acknowledges funds from the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. This work was also made possible by the open-source codes matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), scipy (Jones et al., 2001–), and pandas (McKinney, 2010).