© ESO, 2012

1. Introduction

In spite of their great astrophysical importance (Muzzin et al. 2009), many aspects of asymptotic giant branch (AGB) stars remain poorly understood (e.g. Kerschbaum et al. 2011, and references therein). Among the AGB stars, the luminous Mira variables might be the most enigmatic with their long-period, large-amplitude pulsations and high mass-loss rates. In particular our understanding of their atmospheric structure, which is the link between enriched stellar interior and evolution-dominating wind, is at a turning point, thanks to advances on both the observational side (IR interferometry), and the modelling side with a new generation of self-excited atmosphere models.

The public release of four CODEX model series for M-type Mira variables by Ireland et al. (2011,hereafter ISW11), greatly increases the opportunities for the general community to assess the performance of such dynamic models. Here, we present the first confrontation of the CODEX models (see Sect. 3) to a fully phase-resolved 2 μm spectro-photometric and interferometric dataset covering eight pulsation cycles, as obtained with the Palomar Testbed Interferometer on the Mira TU And. This star has a spectral type M5IIIe, a progenitor mass M = 1.15 − 1.4   M_⊙, metallicity Z = 0.004 − 0.02 (Mennessier et al. 2001), and a GCVS period P ~ 317 d (Samus et al. 2009).

The most intriguing result of this comparison is that there seems to be a systematic discrepancy between the models and our observations around visual maximum: the star appears too hot, and its outer atmosphere too devoid of water. An identification of the source of this discrepancy is hampered by the limited coverage in parameter space of the models, but we discuss possible avenues of further research, and stress the need for larger published model grids and for observational campaigns covering the visual maximum phase in detail.

2. Observations

TU And was extensively observed as part of the large PTI Mira programme (Thompson 2002) during 60 nights from 1999 to 2005, over eight pulsation cycles. The publicly available PTI data archive provides both squared visibility amplitudes (hereafter visibilities) and integrated photon counts. We describe the data reduction and calibration of these observables below.

2.1. Interferometry

The now dismounted Palomar Testbed Interferometer (PTI) was a Y-shaped, long-baseline near-IR interferometer located at Palomar Observatory and operated by JPL/Caltech. It consisted of three 40 cm apertures with separations of 85–110 m and used in pairwise combination (see Colavita et al. 1999, for a detailed description). PTI thus measured only one interferometric observable per spectral channel: the square of the “fringe contrast” or visibility amplitude. Visibilities were recorded in five channels within the K-band (2.0–2.4 μm). Processing of the raw data (correcting for detector biases and read-out noise, and averaging the 125 s scans into fringe amplitudes) was automatically done on site, following procedures described in Colavita (1999). Calibrating for the system visibility, i.e. the response of the measurement system (interferometer + atmosphere) to a point source, was done using the nbCalib tool (Boden et al. 1998) provided by the NExScI, which uses a weighted-mean formalism to estimate the system visibility at the time and position of the science observation from all calibrator measurements that fit the user-defined temporal and spatial constraints (see van Belle et al. 1999, for more details). To ensure calibration stability and consistency, essential for making full use of this data set as an interferometric time series, several calibrators were observed each night and were intercompared by nbCalib.

The stellar parameters of our calibrators were determined by fitting a reddened Kurucz model atmosphere (Kurucz 1993) to high-quality archival photometry. The errors were estimated using Monte Carlo simulations. Table 1 lists the calibrators and their adopted parameters. All calibrators have angular diameters well below the PTI resolution limit of 1 mas (van Belle & van Belle 2005) and are, on average, 10° away from the science target.

The calibrated visibilities of TU And from 50 nights were retained, having a median relative precision ranging from 5% at 2.2 μm, 7.5% at 2.4 μm, and 10% at 2.0 μm.

2.2. Spectro-photometry

Together with the uncalibrated visibilities, PTI delivered the raw photon counts of each observation, typically in blocks of five observations taken within minutes and interleaving science targets with calibrators. We derived the spectral shape over the K-band and flux-calibrated the 2.2 μm channel with a bootstrap method.

The spectral shape calibration procedure consists of dividing the science-target normalised photon counts by the ones of the calibrator and multiplying the result with the intrinsic spectral shape of the calibrator. To obtain an absolute flux in the 2.2 μm channel, one calibrator measurement on either side of the science observation with respect to airmass is required to solve for the atmospheric extinction per unit airmass A_v = (log  [ I₂/I₁ ]  + log  [ S₁/S₂ ] )/(m₂ − m₁) and the instrument spectral response R_v = S₁/(I₁ × exp( − A_vm₁)), with I_i the intrinsic and S_i the measured flux, and m_i the airmass, of the ith calibrator. The calibrated flux is then calculated as I_sci = S_sci/R_v × exp(A_vm_sci), where S_sci denotes the measured photon counts and m_sci the airmass of the science target. For each observing block, the above procedure is applied 5000 times to a randomly picked science observation, corresponding random calibrator observations and calibrator intrinsic fluxes. Finally, the distributions of all observing blocks within a single night are co-added, and the average and standard deviation are taken from the resulting normal distribution, as the final calibrated flux and its error.

To obtain the calibrator intrinsic fluxes, the PTI optical configuration was taken into account: after beam combination and dispersion by a prism, the wavefront was sent through a broadband K Barr filter onto the detector pixels (Colavita et al. 1999). Each of the Kurucz models in the Monte Carlo distributions obtained with the spectral energy distribution fitting, is converted into the effective flux I in each of the PTI channels as (1)with M the Kurucz model, F the filter profile, G a Gaussian of FWHM = 0.7   μm (Colavita, priv. comm.), and λ_l, λ_u the edges of the spectral channel under consideration.

Finally, we constructed a visual light curve by retrieving observed V-band photometry from the AAVSO and AFOEV. Visual phases (hereafter simply phases) were computed by fitting a sine to ten cycles around the times of the PTI observations.

3. The CODEX models

Fig. 1 Spectro-photometry. The two upper panels show the variation in the spectral shape with respect to the 2.2 μm channel, the third panel shows the absolute flux in the latter, and the lower panel the visual flux in magnitudes. In each panel two cycles of the R52 model series are consecutively overplotted as the solid line. On the left, the cycles are shown as a function of visual phase while on the right they are phase-folded.

The CODEX atmosphere models (Ireland et al. 2008) are the most advanced publicly available dynamical models for M-type Mira variables. The pulsations in these models are self-excited and limited in amplitude by the non-linear effects of shocks in the outer atmosphere and by turbulent viscosity. For each model series, a few typical cycles were selected for detailed radiative transfer calculations with an opacity-sampling method in LTE with a phase-resolution of  ~0.1. This results in the centre-to-limb variation (CLV) and spectrum, with Δλ = 0.0002   μm between λ = 0.45−2.50   μm. The current model series are made to match o Cet, R Leo, and R Cas, although even these well-observed Miras have large uncertainties on their basic parameters (e.g. the 2 very different models for R Cas). They were presented to the community and compared to photometric and interferometric observations in ISW11. The model parameters are mass M, luminosity L, metallicity Z, mixing length parameter α_m, and turbulent viscosity parameter α_ν. For a detailed justification of the current set of parameter values, we refer to ISW11.

With a slightly meandering (±5 d) pulsation period of  ~317 d, TU And is compared to the R52 (M = 1.1   M_⊙, L = 5200   L_⊙, Z = 0.02, α_m = 3.5, α_ν = 0.25) and o54 (M = 1.1   M_⊙, L = 5400   L_⊙, Z = 0.02, α_m = 3.5, α_ν = 0.25) model series, which have derived periods of 307 and 330 d, respectively.

When comparing a model series to observations, the only free parameter is the distance, or equivalently, the phase-averaged Rosseland angular diameter θ_R. Considering the angular resolution of PTI, we converted the model intensities to synthetic observables for θ_R ∈  [ 1.5,4.5 ]  mas, in steps of 0.05 mas. Converting the model intensities to the spectro-photometric observables was done with Eq. (1). The squared visibilities are computed following Wittkowski et al. (2004): (2)with F, G, λ_l, and λ_u as in Eq. (1), and (3)the monochromatic Hankel transform of the intensity profile at the baseline B and wavelength λ.

4. Results

Figure 1 shows the resulting photometric time series of TU And as a function of phase, with two consecutive cycles of the R52 model series repeatedly overplotted. The two other available cycles of this series give qualitatively identical results, so are not shown. The Rosseland angular diameter is chosen to match the absolute flux at 2.2 μm, with a resulting value of θ_R = 2.20 ± 0.15 mas or d = 1.13 ± 0.07 kpc. Although the error on the absolute calibration in this channel is rather large, the contamination from upper atmospheric molecular layers at these wavelengths is smallest. This channel thus provides the best view at the continuum photosphere in the spectral region where the bulk of the luminosity is emitted. Moreover, fixing the angular diameter like this gives a consistent picture in both the spectro-photometry and the interferometry.

The R52 model series shows reasonable agreement with the observed visual light curve in the lower panel, considering that (1) the amplitude was not tuned to match this light curve and (2) non-LTE effects play an important role at these wavelengths (Ireland et al. 2008), causing an increase in visual brightness of 0.5 to 1.0 mag. In case the higher value of this correction is valid, the perceived difference between model and observations could point to the presence of an excess in visual light at maximum phases together with a too large model amplitude.

The two upper panels show the variation in the spectral shape with time. A clear discrepancy occurs around times of visual maximum, then the model predicts too much flux. Moreover, the amplitude of this excess shows an interesting wavelength-dependence: the strongest effect is observed at 2.0 μm, a smaller discrepancy at 2.1 and 2.4 μm and a good agreement at 2.3 μm.

Fig. 2 Visibility time series. Each panel shows all squared visibilities as a function of spatial frequency, obtained at the visual phase given in the lower left corner of the panel. Different symbols are used for observations from different cycles. The same two model cycles as in Fig. 1 are overplotted as the full and dashed lines in each panel.

The visibility time series is presented in Fig. 2, overplotted with the same cycles of the R52 model series. Each panel displays the visibilities at a given phase, using different symbols for observations coming from different cycles (at least two per panel). The model reproduces the observations well at phases 0.3–0.6, but not around visual maximum. Comparison with other cycles does not give better results. In fact, this reveals that the model has larger cycle-to-cycle variations at pre-maximum phases than are seen in the observations. A closer look at the wavelength dependence shows that the discrepancy is least pronounced at 2.3 and 2.4 μm. Between 2.0 and 2.2 μm, the model shows very little wavelength-dependence, unlike the observations which show a clear decrease in visibility towards 2.0 μm.

This observation is supported by the more quantitative assessment presented in Figs. 3 and 4. Epochs with observations on at least two different baselines were selected, to resolve the degeneracy between CLV basic size and shape. Each of these data sets with good uv-coverage was compared to each model, independently of model phase and for a range of θ_R. The figures show the resulting χ²-maps. In all panels a curved valley shows up: every data set can be fitted with a model of any phase, albeit with a significantly larger diameter for visual maximum models. However, within these valleys a better χ² is typically found with a model at the correct (photometric) diameter but with a phase nearer to visual minimum than the actual value of the observations. In conclusion, we find that the models around visual maximum should resemble those of visual minimum more in CLV shape and wavelength-dependence.

Finally, we also compared the observations to the similar o54 model series, except for its luminosity. The results for o54 are alike and thus not shown here. The visibilities seem to match slightly better, but the light curve discrepancies become larger.

5. Discussion

Fig. 3 The χ2 maps of two data sets with a good uv-coverage (see text). The log χ2 is shown as a function of Rosseland angular diameter and phase, of the same two model cycles as in Figs. 1 and 2. The phase of the data set is given in the upper right corner of the panel and is plotted as the horizontal dotted line. The vertical dotted line marks the photometric diameter. The intersection of the dotted lines shows where the minimum χ2 should be if the model is correct.

Ireland et al. (2008) demonstrate that the CODEX models have reached a level of maturity, concerning computational issues and input physics, that allows quantitative comparisons with observations. Pulsations are self-excited and molecular layers are produced naturally and in a physically self-consistent way.

Nevertheless, our observations of TU And show a clear discrepancy with the models around phases of visual maximum, both in the spectrum and in the visibilities. The most natural interpretation of the wavelength-dependence of this discrepancy seems to be in terms of the shape of the H₂O opacity curve, which shows an increase from 2.2 μm towards 2.0 and 2.4 μm. The formation of H₂O layers depends sensitively on the density, temperature, and pressure stratification in the upper atmosphere. In the models, the conditions in the upper atmosphere are set by the dynamics. Shock fronts that emerge from the continuum photosphere between phases –0.3 and –0.1 and then propagate outwards lead to the formation of an H₂O layer around phase 0.2–0.3, which then strongly decreases in strength at phase 0.7–0.8. At φ_V = 0.5 both visibilities and spectra show good agreement. However, as demonstrated by Figs. 1–3, the near-maximum observations are represented better by model phases nearer to minimum, and thus with the presence of an H₂O layer in the atmosphere. This interpretation is supported by the better agreement at 2.3 and 2.4 μm. There CO is an important source of opacity, in the photosphere and in the extended atmosphere, even at visual maximum. As a much more stable molecule, CO depends less on the exact conditions and partly hides the lack of an H₂O molecular layer in the spectrum and visibilities.

In retrospect, it is interesting to note that the image reconstruction of T Lep (P ~ 380 d, Le Bouquin et al. 2009), made at a phase of 0.8–0.9, shows the clear signature of an extended molecular layer with a brightness ratio of  ~10%, unlike its closest model series o54, which supports our findings. Woodruff et al. (2009) found fair agreement between the wavelength dependence of their measured 1.1–3.8 μm UDs and the predictions of the CODEX models, for several Miras at phases 0.5 to 0.8, but needed slightly larger model continuum diameters than predicted. On the other hand, Wittkowski et al. (2011) have recently compared AMBER medium spectral resolution observations of four Miras to the CODEX models and found good agreement, although they compared the models to stars of a different period, hence luminosity.

In summary, at visual maximum the current models predict a continuum photosphere that is too compact and hot as already suggested by ISW11, in combination with too small an H₂O optical depth in the outer atmosphere, as found here.

The basic question now is, can this be salvaged within the current modelling framework? Is it just a question of tuning the parameters to more realistic values? Following the discussion of ISW11, it seems difficult to find any combination of parameter values that fits all observables. Accepting the LMC period-luminosity relation and assuming that the effect of metallicity is insufficient (see ISW11), the only free parameters are M, α_m, and α_ν. Since the last mainly influences the pulsation amplitude, but not the period or the “parent star” radius, the above authors conclude that the problem at hand is one of calibrating the mixing length parameter as a function of mass. Decreasing the

rather high α_m = 3.5 to a more “realistic” value requires an increase in mass or a decrease in metallicity (albeit with a smaller dependency) to preserve the period, which would in case of o Cet be inconsistent with its kinematically determined age. TU And though, was classified as a member of a population with initial masses in the 1.15 − 1.4   M_⊙ range and Z ∈  [ 0.004,0.02 ]  (Mennessier et al. 2001), which makes it a better candidate to test the above hypothesis.

We can at this point not gauge the impact of using an α_ν which is not tuned to the amplitude of TU And. As a result, we feel that the model grid should span some scope in this parameter to allow more comparisons with “small-amplitude” pulsators.

Could either change in the model parameters also induce a stabilizing effect on the H₂O layer, as required by our observations, or is there a need to add more physics to the models?

6. Conclusions

The combination of a spectro-photometric and interferometric time series is a powerful tool for testing the validity of pulsating AGB model atmospheres, even at low spectral resolution. PTI observed a large sample of M-type Miras and semi-regulars, several with a similar sampling to the one presented here, covering periods between 150 and 470 d. To fully exploit this rich database, better coverage in the model parameter space is urgently needed, even though it requires a large computational effort. Only in this way can the model parameters be calibrated and (simply) linked to Mira observable properties. From the observational point of view, interferometric studies with better uv-coverage and spectral resolution should focus on phases 0.0 to 0.15, since none of the current model series predicts any extended molecular atmosphere in that phase range.

Online material

Fig. 4 Similar to Fig. 3, but for the four other data sets with good uv-coverage. The log χ2 is shown as a function of Rosseland angular diameter and phase, of the same two model cycles as in Figs. 1 and 2. The phase of the data set is given in the upper right corner of the panel and is plotted as the horizontal dotted line. The vertical dotted line marks the photometric diameter. The intersection of the dotted lines shows where the minimum χ2 should be if the model is correct.

Acknowledgments

The Palomar Testbed Interferometer was operated by the NASA Exoplanet Science Institute and the PTI collaboration. It was developed by the Jet Propoulsion Laboratory, California Institute of Technology with funding provided from the National Aeronautics and Space Administration. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research.

References

All Tables

All Figures

Fig. 1 Spectro-photometry. The two upper panels show the variation in the spectral shape with respect to the 2.2 μm channel, the third panel shows the absolute flux in the latter, and the lower panel the visual flux in magnitudes. In each panel two cycles of the R52 model series are consecutively overplotted as the solid line. On the left, the cycles are shown as a function of visual phase while on the right they are phase-folded.

In the text

	Fig. 2 Visibility time series. Each panel shows all squared visibilities as a function of spatial frequency, obtained at the visual phase given in the lower left corner of the panel. Different symbols are used for observations from different cycles. The same two model cycles as in Fig. 1 are overplotted as the full and dashed lines in each panel.
In the text

Fig. 3 The χ2 maps of two data sets with a good uv-coverage (see text). The log χ2 is shown as a function of Rosseland angular diameter and phase, of the same two model cycles as in Figs. 1 and 2. The phase of the data set is given in the upper right corner of the panel and is plotted as the horizontal dotted line. The vertical dotted line marks the photometric diameter. The intersection of the dotted lines shows where the minimum χ2 should be if the model is correct.

In the text

Fig. 4 Similar to Fig. 3, but for the four other data sets with good uv-coverage. The log χ2 is shown as a function of Rosseland angular diameter and phase, of the same two model cycles as in Figs. 1 and 2. The phase of the data set is given in the upper right corner of the panel and is plotted as the horizontal dotted line. The vertical dotted line marks the photometric diameter. The intersection of the dotted lines shows where the minimum χ2 should be if the model is correct.

In the text