© V. Vasilyev et al. 2019

1 Introduction

Oxygen is the third most abundant chemical element and the most abundant metal in the Universe. As with other α-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; 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 α-, 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, and 845 nm triplets can be usedto 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 and 636.4 nm lines, or to the permitted high-excitation O I 615.8 nm triplet. However, the near-IR triplets are sensitive to departures from LTE. For a given oxygen abundance, taking into account non-LTE effects tendsto 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 differentspatial 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 κ-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 Sect. 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 nm triplet. The results including line profiles, abundance errors and biases of the standard approach are presented in Sect. 3. In Sect. 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.

Fig. 1 Light curves in terms of the emergent bolometric flux of models M3 (left panel) and M9 (right panel) which exhibit pulsational periods of 3 and 9 days, respectively. Blue dots indicate the instances in time for which 2D non-LTE spectral syntheses were performed.

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_eff = 5600 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 and 6730 K; in particular, during the seventh cycle the maximum effective temperature reaches a value above 6700 K.

While the 2D models have a depth-independent gravitational acceleration g_o, spectral lines are in fact sensitive to the effective gravitational acceleration, $g_{\text{eff}}=g_{\text{o}}-\frac{\text{d}v}{\text{d}t}$, 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 × N_z = 600 × 500 along the horizontal and vertical directions. The geometrical extensions of the models along the horizontal and vertical directions are l_x × l_z = 1.1 × 10¹² × 0.5 × 10¹² cm and l_x × l_z = 3.6 × 10¹² × 1.7 × 10¹² 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_{\text{eff}}^{\text{2D}}$ was calculated from the bolometric radiative flux passing the top boundary at each instance in time. The kinematic acceleration $\frac{\text{d}v}{\text{d}t}$, which is variable during the pulsational cycle, leads to a variation in g_eff in the range ≈(0.6…3.5)g_o for M9, and ≈(0.5…3.0)g_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_eff and log g_eff over the pulsational cycles for models M3 and M9 in Fig. 2.

Fig. 2 Grid of 1D hydrostatic models (red dots) in the log g–Teff plane. The 2D snapshots defined individually by log geff and Teff for model M3 (left panel) and M9 (right panel) are indicated by blue dots.

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 α was fixed to 1.5 for all 1D models. In order to cover the parameters encountered in M3 and M9, T_eff, and log g for the 1D models were chosen in the range 5200 < T_eff∕K < 6800 with a step ΔT_eff = 100 K and 1.1 < log g < 2.9 with Δ log g = 0.2, respectively. We depict the grid of 1D models in the log g–T_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 calculate2D 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 Δτ_R ≲ 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; Kovtyukh et al. 2005, and the later studies of Vasilyev et al. 2017; Proxauf et al. 2018).

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(v_2D) = W(ξ_t), where W(v_2D) is calculated with the full 2D velocity field v_2D, and W(ξ_t) by replacing v_2D by a depth-independent isotropic microturbulent velocity ξ_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 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 ξ_t was introduced, and was varied between 0.0 and 4.0 km s⁻¹ with a step 1.0 km s⁻¹. Calculations were performed for several oxygen abundances − 0.8 < [O∕Fe] < +0.8 with a step 0.4 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, $${\Delta }_{\text{LTE}}^{\text{NLTE}}\left(\log {\varepsilon }^{\text{LTE}}\right)=\log {\varepsilon }^{\text{NLTE}}-\log {\varepsilon }^{\text{LTE}},$$(1)

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 }_{\text{1D}\text{NLTE}}^{\text{2D}\text{NLTE}}$ and ${\Delta }_{\text{1D}\text{LTE}}^{\text{2D}\text{NLTE}}$, we can look at the variation of the departure coefficients, b(τ_R) = n_NLTE∕n_LTE, for the 1D hydrostatic models and 2D models. For the first, b^1D values were calculated for the lower 3s⁵ S and upper 3p⁵ P levels of the O I 777 nm triplet in the 4D parameter space τ_R, T_eff, log g_eff, ξ_t. For each combination ${[T_{\text{eff}},\log g_{\text{eff}},{\xi }_{\text{t}}]}^{\text{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(τ_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^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 τ_R and photometricphase for the lower and upper levels are shown by solid blue and red lines, respectively. One-dimensional interpolated b^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 nm triplet, namely 3s⁵ S and 3p⁵ P, are almost always overpopulated during the pulsating cycle. At photometric phases ϕ_ph ≈ 0.3…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 }_{\text{1D}\text{NLTE}}^{\text{2D}\text{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 nm triplet profiles

We illustrate the temporal evolution (over a pulsation cycle) of the normalised emergent LTE and non-LTE O I 777 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 ϕ_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_eff spans up to around 1000 K, corresponding to around a 0.3–0.5 dex change in line strength; instead, at fixed T_eff, log g spans up to around 0.4 dex, corresponding to a less than 0.04 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.

Fig. 3 1D b1D(τR) and median 2D b2D(τR) departures as a function of τR for the lower 3s5 S (blue) and upper 3p5 P (red) levels (dashed and solid lines, respectively) for a particular cycle. The median 2D departures are results of the horizontal average along the surfaces τR = const. Departures from LTE b1D(τR) are based on interpolation of results of calculations using 1D hydrostatic models by fixing ${[T_{\text{eff}},\log g_{\text{eff}},{\xi }_{\text{t}}]}^{\text{2D}}$.

Fig. 4 Variation in the O I 777 nm triplet during a particular pulsational cycle for LTE (blue) and non-LTE (red) cases calculated for the three-day M3 Cepheid model (left panel) and for the nine-day M9 Cepheid model (right panel). The dashed vertical lines indicate the lines’ rest wavelengths.

Fig. 5 Equivalent widths for the O I 777 nm triplet integrated over all components as functions of effective temperature and effective surface gravity for [O∕Fe] = −0.4, 0.0, +0.4 and models M3 and M9.

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 temperatureand 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 ratiosin 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_eff ≈ 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 nm component of the O I 777 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 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_eff − log g −Amplitude −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 g_eff of an Lagrangian mass element located near optical depth unity in the 2D models. This differs from a g_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.

Fig. 6 2D NLTE vs. 2D LTE abundance corrections for the O I 777 nm triplet integrated over all components as functions of effective temperature and effective surface gravity for [O∕Fe] = −0.4, 0.0, +0.4 and models M3 and M9.

Fig. 7 1D NLTE vs. 1D LTE abundance corrections for the 777.4 nm component of the O I 777 nm triplet as a function of effective temperature. Literature results of Luck et al. (2013) based on 1D ATLAS9 model atmospheres are shown as black stars. Our 1D results are shown as dots taking effective temperatures and gravities from the three-day (red) and nine-day (blue) Cepheid model.

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 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_{\text{eff}},\log g_{\text{eff}},{\xi }_{\text{t}}\right]}^{\text{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 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 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, wedo 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_eff ≲ 5700 K, the abundance corrections are almost insensitive to the change in T_eff and are around − 0.3 dex, ranging from − 0.2 to − 0.4 dex. For higher effective temperatures the abundance corrections become more severe, decreasing to − 1.5 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, g_eff drops by a factor of three while T_eff increases by 130 K. For this short time interval, the increase in the O I 777 nm triplet component equivalent widths in LTE and non-LTE cases are primarily due to the sensitivity to the changes in log g_eff. In our calculations, the relative changes of LTE and non-LTE equivalent widths are ${\left(\frac{\Delta W}{W}\right)}_{\text{LTE}}\approx 0.4$ and ${\left(\frac{\Delta W}{W}\right)}_{\text{NLTE}}\approx 0.2$, respectively. As a result, the 2D non-LTE versus 1D LTE abundance error increases from − 1.35 to − 1.1 dex.

Fig. 8 2D non-LTE vs. 1D non-LTE (upper panels) and 2D non-LTE vs. 1D LTE (lower panels) abundance corrections for the O I 777 nm triplet integrated over all components as functions of effective temperature for the nine-day (left panels) and three-day (right panels) Cepheid model. The photometric phase ϕph is indicated by colour. Solid and dashed black lines depict the temperature dependence for M9 and M3, respectively, when fittingthe numerical data by a parabola.

4 Discussion

4.1 Optimal phases for 1D analyses

In practice, spectroscopic analyses of the O I 777 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 ± 0.1 dex, and in the case of 1D non-LTE analysis 0 ± 0.1 dex at photometric phases ϕ_ph ≈ 0.3…0.8 independentof 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 phaseswould 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 < ϕ_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σ uncertainties. For the data obtained at all phases the slope is − 0.053 ± 0.003 dex kpc⁻¹, whereas with data at optimal phases we get a slightly steeper gradient − 0.056 ± 0.004 dex kpc⁻¹. 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 metallicities, 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 Minchev et al. (2013) 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 O I 777 nm triplet owing to its temperature sensitivity, which is amplified by temperature-sensitive 3D non-LTE effects. We expect the dependencesto be stronger in Cepheids with higher mean effective temperatures, to be more metal-rich, and to have longer periods.

Fig. 9 Dependence of the 2D non-LTE vs. 1D non-LTE (upper panel) and 2D non-LTE vs. 1D LTE (lower panel) abundance corrections for the O I 777 nm triplet integrated over all components as a function of photometric phase for the model M3 (red) and M9 (blue).

Fig. 10 Plot of 1D non-LTE oxygen abundances in the Galactic disc, as presented by Korotin et al. (2014). Abundances that are based on the analysis of spectra taken at all phases and at optimal phases 0.3 < ϕph < 0.8 are shown as light blue dots and red circles, respectively. Linear fits for these two cases including 1σ uncertainties are presented in red and blue lines.

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 }_{\text{2D}\text{LTE}}^{\text{2D}\text{NLTE}}$ 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 ϕ_ph ≈ 0.3…0.8, independentof 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_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

V.V., H.G.L., and B.L. 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. A.M.A. 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).

References

All Figures

	Fig. 1 Light curves in terms of the emergent bolometric flux of models M3 (left panel) and M9 (right panel) which exhibit pulsational periods of 3 and 9 days, respectively. Blue dots indicate the instances in time for which 2D non-LTE spectral syntheses were performed.
In the text

	Fig. 2 Grid of 1D hydrostatic models (red dots) in the log g–Teff plane. The 2D snapshots defined individually by log geff and Teff for model M3 (left panel) and M9 (right panel) are indicated by blue dots.
In the text

Fig. 3 1D b1D(τR) and median 2D b2D(τR) departures as a function of τR for the lower 3s5 S (blue) and upper 3p5 P (red) levels (dashed and solid lines, respectively) for a particular cycle. The median 2D departures are results of the horizontal average along the surfaces τR = const. Departures from LTE b1D(τR) are based on interpolation of results of calculations using 1D hydrostatic models by fixing ${[T_{\text{eff}},\log g_{\text{eff}},{\xi }_{\text{t}}]}^{\text{2D}}$.

In the text

	Fig. 4 Variation in the O I 777 nm triplet during a particular pulsational cycle for LTE (blue) and non-LTE (red) cases calculated for the three-day M3 Cepheid model (left panel) and for the nine-day M9 Cepheid model (right panel). The dashed vertical lines indicate the lines’ rest wavelengths.
In the text

	Fig. 5 Equivalent widths for the O I 777 nm triplet integrated over all components as functions of effective temperature and effective surface gravity for [O∕Fe] = −0.4, 0.0, +0.4 and models M3 and M9.
In the text

	Fig. 6 2D NLTE vs. 2D LTE abundance corrections for the O I 777 nm triplet integrated over all components as functions of effective temperature and effective surface gravity for [O∕Fe] = −0.4, 0.0, +0.4 and models M3 and M9.
In the text

	Fig. 7 1D NLTE vs. 1D LTE abundance corrections for the 777.4 nm component of the O I 777 nm triplet as a function of effective temperature. Literature results of Luck et al. (2013) based on 1D ATLAS9 model atmospheres are shown as black stars. Our 1D results are shown as dots taking effective temperatures and gravities from the three-day (red) and nine-day (blue) Cepheid model.
In the text

Fig. 8 2D non-LTE vs. 1D non-LTE (upper panels) and 2D non-LTE vs. 1D LTE (lower panels) abundance corrections for the O I 777 nm triplet integrated over all components as functions of effective temperature for the nine-day (left panels) and three-day (right panels) Cepheid model. The photometric phase ϕph is indicated by colour. Solid and dashed black lines depict the temperature dependence for M9 and M3, respectively, when fittingthe numerical data by a parabola.

In the text

	Fig. 9 Dependence of the 2D non-LTE vs. 1D non-LTE (upper panel) and 2D non-LTE vs. 1D LTE (lower panel) abundance corrections for the O I 777 nm triplet integrated over all components as a function of photometric phase for the model M3 (red) and M9 (blue).
In the text

	Fig. 10 Plot of 1D non-LTE oxygen abundances in the Galactic disc, as presented by Korotin et al. (2014). Abundances that are based on the analysis of spectra taken at all phases and at optimal phases 0.3 < ϕph < 0.8 are shown as light blue dots and red circles, respectively. Linear fits for these two cases including 1σ uncertainties are presented in red and blue lines.
In the text