© ESO 2018

1 Introduction

Most of the chemical elements in the Universe were forged inside evolved stars. As one of the latest stages in the evolution of massive stars, red supergiant (RSG) stars contribute to this enrichment through their mass loss. The mechanism that launches the material away from the star remains unknown. One scenario involves convection. From spectroscopic observations, Josselin & Plez (2007) suggested that by lowering the effective gravity, the turbulent velocity field associated with convection allows the radiative pressure on molecular lines to start the outflow. Schwarzschild (1975) predicted that photospheric convection on RSG stars would be different from what is known on solar-type stars: only a handful of giant granules would be present on the stellar surface. Near-infrared (NIR) observations by Haubois et al. (2009) on Betelgeuse, the prototypical M-type RSG star, showed that only one or two bright and large hot spots were present on the photosphere. These features were interpreted as the top of convective granules. Additional observations were recently obtained that point toward a convection-based launch mechanism: ALMA observations show a large bright spot on the photosphere of Betelgeuse (O’Gorman et al. 2017). Its position matches the direction of a strong linearly polarized clump observed with VLT/SPHERE at three stellar radii (Kervella et al. 2016) that has been interpreted as recently formed dust. Kervella et al. (2018) determined that the rotation axis of the star is also aligned with these two features and suggested that enhanced mass loss was emitted from the polar region due to the long term presence of a “rogue” convective cell. However, the velocity field derived by Ohnaka et al. (2017) on the RSG Antares leads them to the conclusion that convection alone is not able to explain the atmospheric extension and motions of this star. A similar conclusion was obtained by Arroyo-Torres et al. (2015) on a sample of RSG stars.

Furthermore, Chiavassa et al. (2011b) showed that these large granules can additionally cause photocenter displacements significant enough to bias parallax measurements. As RSG stars can be used as bright candles at large distances, this effect has important consequences.

Montargès et al. (2016) monitored the H-band photosphere of Betelgeuse between January 2012 and November 2014. They obtained four epochs of observations showing a large feature (characteristic size ~R_⋆) departing from spherically symmetric limb-darkened disk (LDD) models. These observations were correlated with January and November 2014 spectropolarimetric measurements at optical wavelengths obtained at the Narval instrument mounted on the Telescope Bernard Lyot (TBL) at the Pic du Midi observatory. From these observations, Aurière et al. (2016) discovered the linearly polarized spectrum of Betelgeuse. Their analysis points to a continuum depolarization of Betelgeuse (due to scattering at photospheric level) that may be related to brightness inhomogeneities lying at the stellar surface. The location and evolution of these inhomogeneities were mapped with an analytic model.

We report here NIR interferometric observations conducted on the RSG star CE Tau (119 Tau, HR 1845, HD 36389). CE Tau has a M2Iab-Ib spectral type. It has no reported companion. Cruzalèbes et al. (2015) observed significant departure from centrosymmetry on this star using the AMBER instrument at the Very Large Telescope Interferometer (VLTI). In Sect. 2 we present the observations we obtained with VLTI/PIONIER. The angular diameter of the star is determined in Sect. 3. Photospheric features are studied in Sect. 4 using several approaches: classical spotty models, image reconstruction, and 3D radiative hydrodynamics simulations. The updated angular diameter value is used to derive updated stellar parameters. The contrast of the convective pattern is discussed in Sect. 5. We present a summary and conclusions in Sect. 6.

Fig. 1 (u, v) coverage of our three epochs of VLTI/PIONIER observations of CE Tau. Left: 14 November 2016. Center: 22 November 2016. Right: 23 December 2016. The underlying image corresponds to the best fit LDD model of Sect. 3

2 Observations and data reduction

CE Tau was observed on 14 and 22 November and on 22 and 23 December 2016 at the European Southern Observatory’s VLTI (Haguenauer et al. 2010) located on top of Cerro Paranal in Chile. We used the Precision Integrated-Optics Near-infrared Imaging ExpeRiment (PIONIER; Le Bouquin et al. 2011) instrument, equipped with the RAPID detector, fed by the four 1.8 m diameter Auxiliary Telescopes (AT) in their compact (stations A0-B2-C1-D0) and intermediate (stations D0-G2-J3-K0) configurations. Ground baselines were between 6.9 and 96.9 m. The GRISM was set in the optical path of PIONIER, providing its highest available spectral resolution of R ~ 30, delivering six spectral channels between 1.51 and 1.77 μm. The log of the observations is available in Table A.1. We note that other executions of the observing blocks exist in the archive, but they had a low quality grade. The data from 22 December show some visibility loss. As the observation was repeated on 23 December without any issue, and because we do not expect the photosphere of a RSG star to evolve significantly over 24 h, we will not consider the 22 December data further. The (u, v) coverage for each of the three epochs we consider is represented in Fig. 1.

The data were reduced and calibrated with the publicly available PIONIER pipeline (Le Bouquin et al. 2011) version 3.51. The angular diameters of our interferometric calibrators are listed in Table 1. The uncertainties on the observables are directly computed by the pipeline: on the uncalibrated data it derives the statistical dispersion over 100 scans of each ~30 s exposure. Then for the calibrated product it quadratically adds the error from the transfer function. Each individual observation results in six squared visibilities and four closure phases per spectral channel.

3 Angular diameter measurements

Our data from 14 November 2016 were acquired in the compact configuration. These observations cover only the first lobe of the visibility function (Fig. 1, left). Data from 22 November and 23 December cover the first three lobes (Fig. 1, center and right). When plotting the visibility versus the spatial frequency (Fig. 2), we observe that the low spatial frequency first lobe data between 22 November and 23 December are invariant (while the signal changes for the longer baselines). Lower spatial frequencies probe larger scale features andin the case of the first lobe, the general stellar shape. These are crucial for reliable model fitting and image reconstruction. Therefore, we decided to merge the 14 November data on the compact configuration (very short baselines) with the intermediate configuration data of 22 November (November dataset hereafter) and with the 23 December data (December dataset hereafter). This allowed us to have short baseline data on both epochs.

To determine the angular diameter of the star, we used both a uniform disk (UD, I = I₀) and a LDD power-law model (I∕I₀ = μ^α). The visibility amplitude of the latter is given by Hestroffer (1997) $$V_{\text{LDD}}(s)=\text{Γ}(\nu +1)\frac{J_{\nu }(x)}{{(x/2)}^{\nu }},$$(1)

where ν = α∕2 + 1, s is the spatial frequency, x = πsθ_LDD, θ_LDD the LDD angular diameter, J_ν is the first species Bessel function of order ν, and Γ is the Euler function.

In order to derive the angular diameter, we consider here only the 1.62 μm channel, closest to the H⁻ opacity minimum (Gray 2008). A global fit showed important deviations in the third lobe, probably due to contaminations by inhomogeneities. To avoid a bias in the diameter estimation, we present only the fit restricted to the first and second lobes (spatial frequencies lower than 220 arcsec⁻¹). The results of the fit are given in Table 2 and the LDD model is represented in Fig. 2 with the data. The UD model reproduces the data poorly and will not be considered further.

The fits of the two epochs of observations by the LDD model give different but compatible values. For a fair comparison with the models including inhomogeneities presented in Sect. 4, we also derived the reduced χ² for the best LDD model taking into account both the visibilities and the closure phases of the entire spatial frequency domain. We obtain ${\tilde{\chi }}^2=3106$ and 4578 for the November and December datasets, respectively. We do not fit the LDD model to the closure phases as this observable constrains disk models poorly (our dataset does not sample the flips from 0 to π exactly). The most recent angular diameter measurement of this star is given by Cruzalèbes et al. (2013) with 9.97 mas in the K band measured with VLTI/AMBER on November 2009. These authors compiled all the previous angular diameter measurements of this star. For different spectral domains the values range between 9.4 and 13.0 mas for indirect methods, 9.1 and 17 mas for lunar occultations, and 9.3 and 10.68 mas for long baseline interferometry. Therefore, with values of 10.09 ± 0.09 and 10.18 ± 0.07 mas, our LDD diameters appear compatible with these results. The poor reduced χ² of the LDD model when compared to the closure phase, although the squared visibility value is acceptable, suggests the presence of photospheric structures that cannot be reproduced by simple disk models.

Fig. 2 Squared visibility in the first two lobes measured by VLTI/PIONIER on CE Tau at 1.62 μm. The purple squares correspond to the November dataset. The orange dots correspond to the December dataset. The continuous curves are the best fit LDD power-law models for both epochs.

4 Photospheric features

4.1 Limb-darkened disk and Gaussian spots

4.1.1 Model with one spot

The presence of inhomogeneities on the photosphere of RSG stars is expected. Previous observations (e.g., Haubois et al. 2009; Chiavassa et al. 2010; Baron et al. 2014; Montargès et al. 2016; Aurière et al. 2016; Ohnaka et al. 2017; Wittkowski et al. 2017) have shown that they are common. They are also predicted by models (Schwarzschild 1975; Chiavassa et al. 2011a). The most economical models for inhomogeneities are spots, either uniform disk or Gaussian. Here, following the models used by Montargès et al. (2016); Montargès et al. (2017), we use Gaussian spots. We consider the LDD model presented in Sect. 3, and we add a Gaussian spot at the position (x_center, y_center) relative to the center of the stellar disk. We denote w_LDD and w_spot the peak flux of the LDD model and the spot, respectively, and FWHM is its full width at half maximum. We normalize the model using $$w_{\text{LDD}}+w_{\text{spot}}=1.$$(2)

The complex visibility of the model is then $$V_{\text{model}}=w_{\text{LDD}}{V}_{\text{LDD}}+w_{\text{spot}}{V}_{\text{spot}},$$(3)

with $$\begin{array}{lll}\hfill & V_{\text{spot}}(u,v)\hfill & =\exp \left[-\frac{{(2\pi f\sigma )}^2}{2}\right]\hfill \\ \hfill & \hfill & \times \exp \left[-2i\pi (u{x}_{\text{center}}+v{y}_{\text{center}})\right],\hfill \end{array}$$(4)

where i² = −1, $f=\sqrt{u^2+v^2}$ and $\sigma =FWHM/(2\sqrt{2\ln (2)})$.

Modeling such features is difficult as the χ² distribution becomes very complex, as discussed by Baron et al. (2014). Therefore, the fitting process is done in two steps for each epoch. In the first step, the parameters of the limb-darkened disk are fixed to the values derived in Sect. 3. The FWHM and w_spot explore the ranges [0.05; 5 mas] by steps of 0.5 mas and [−0.5; 0.5] by steps of 0.05, respectively. For each couple (FWHM, w_spot) a χ² map is built by deriving the reduced χ² associated with the model for various positions of the Gaussian spot on the stellar disk. We used 50 × 50 pixels maps whose edges correspond to the stellar radius. The maps obtained for the best couple (FWHM, w_spot) are shown in Fig. B.1. The parameters associated with the minimum χ² over this set of maps are used as initial guesses in the second step of the fitting process. In this next step, a Levenberg–Marquardt fit is done on all parameters. The results are presented in Table 3. We also derived the F parameter: $$F=\frac{{\chi }_{\text{LDD}}^2-{\chi }_{\text{1spot}}^2}{N_{\text{param},\text{1spot}}-N_{\text{param},\text{LDD}}}\times \frac{N_{\text{data}}-N_{\text{param},\text{1spot}}}{{\chi }_{\text{1spot}}^2}.$$(5)

This parameter allows us to determine the significance of a fit with a large number of parameters. For a (2, 6) distribution and a significance of 5%, the F parameter must be above 5.143¹. In principle, the F parameter should be derived separately for the squared visibilities and the closure phases, but Montargès et al. (2017) have shown that it is irrelevant for the bright targets accessible to optical interferometers.

The very high values of the F parameter indicate that adding a single Gaussian spot strongly improves the fitting process without overfitting. We note that with our approach, bright and dark spots were allowed (positive or negative flux), but the model converged to a bright spot, with very nearby parameters for both epochs.

Attempts were made to fit a two-spot model, but no convergence was reached for the December dataset. Our conclusion is that the χ² distribution for a two-spot model is much more complex: there is an important correlation between the characteristics of the two spots, making the determination of the absolute maximum an extremely complex task. Therefore, we limit our analysis to the single-spot case.

4.1.2 Basic stellar parameters

From the determination of the LDD angular diameter, we are able to derive several basic characteristics of the star. With the parallax of 1.82 ± 0.26 mas (van Leeuwen 2007), we are able to derive the linear radius of the star R = 587 ± 85 R_⊙ in November, and R = 593 ± 86 R_⊙ in December. These values are in agreement with the value of 601 ± 83 R_⊙ from Cruzalèbes et al. (2013).

To determine the bolometric flux, we can use photometric measurements from Lee (1970) and Ducati (2002) ranging from the U to the N band (0.36–10.2 μm). After determining the difference between the B-V of these data (2.08) and the intrinsic B–V color of a M2Iab-Ib star (1.69) reported in Elias et al. (1985), we can compute the interstellar extinction A_λ in each filter using the data in Savage & Mathis (1979). From this we derive the bolometric flux F_UBVRIJHKLN = 7.01 × 10⁻⁹ W m⁻². Because CE Tau is a semi-regular variable, we use the 1 magnitude visual dispersion reported by the American Association of Variable Star Observers (AAVSO) to adopt a 12.8% uncertainty on the bolometric flux (Ohnaka et al. 2013): ± 8.98 × 10⁻¹⁰ W m⁻².

The LDD angular diameter and the bolometric flux give us access to the effective temperature: T_eff = 3820 ± 135 K in November and T_eff = 3801 ± 134 K in December. The 3σ range of these values largely encompasses the value of 3660 K given for a M2 RSG by Levesque et al. (2005) and the 3700 K of Luck & Bond (1980). We can also use the parallax and the bolometric flux to derive the luminosity of the star: $\log L/L_{\odot }={4.82}_{-0.16}^{+0.12}$, a value compatible with that found by Cruzalèbes et al. (2013), who derived a value of 4.63 with an uncertainty of 13%.

Placing CE Tau in a Hertzsprung–Russell (H–R) diagram (Fig. 3), we see that the derived values agree remarkably well with the evolutionary track (Ekström et al. 2012) of a 15 M_⊙ star with rotation (the rotation rate on the zero age main sequence is v_ini∕v_crit = 0.4) and solar metallicity (Z = 0.014). According to Luck & Bond (1980) the metallicity of CE Tau is [Fe/H] = 0.11. Given the error bars on the parameters we adopt an uncertainty of 2 M_⊙ on the stellar mass. For an initial mass of 15 M_⊙, the evolutionary models predict that at its current age, CE Tau weighs $M_{\text{cur}}={14.37}_{-2.77}^{+2.00}$ M_⊙. With this mass and the linear radius, we can derive the surface gravity: $\log g={0.05}_{-0.17}^{+0.11}$, a value compatible with logg = 0.07 obtained by Luck & Bond (1980). We can also use the location of CE Tau on the H-R diagram and the evolutionary models (with and without rotation) to estimate the age of the star: ${13.9}_{-2.5}^{+1.0}$ Myr.

Fig. 3 Hertzsprung–Russell diagram with the position of CE Tau marked by the red cross. The continuous (respectively, dashed) lines correspond to the non-rotational (respectively, rotational) evolutionary tracks of Ekström et al. (2012) for solar metallicity abundances (Z = 0.014).

4.2 Image reconstructions

As fitting analytical models to reproduce photospheric inhomogeneities does not give unequivocal results due to the complexity of the χ² distribution, we decided to continue our analysis with the image reconstruction technique. Our (u, v) coverage actually allows this procedure if we take into account the synthesized beam, while our limited resolution relative to the angular diameter of the star does not enable us to use the statistical approach of Montargès et al. (2017).

To produce images from the datasets, we use SQUEEZE, a compressed sensing-based image reconstruction tool. SQUEEZE uses a Markov chain Monte Carlo (MCMC) approach to search the imaging probability space via parallel tempering (Baron et al. 2010). This method allows the simultaneous use of a variety of regularizers, including the l0 norm non-convex regularizer.

Because of the limited (u, v) coverage of our observations, we decided to use all the spectral channels at once and to start our reconstruction process by determining the best regularizers and regularization strengths for our data. We did this by using the 3D radiative hydrodynamics (RHD) model described in Sect. 4.3 as the source image for a simulated observation, as this allowed us to compare reconstruction results to a known source (center and bottom images of Fig. 5). We produce the simulated observation using OIFITS-SIM², copying the (u, v) coverage and noise statistics of our December observations. The reconstructions use masks and initial images based on uniform disks of sizes 10–12 mas, depending on the pixel scale. In order to correct for artifacts in the image reconstruction process, we ran five MCMC chains, each consisting of 500 iterations, in the end producing an average and error image after coalligning the mean result of each chain using the subpixel registration algorithm of Guizar-Sicairos et al. (2008). In order to test for possible super-resolution in the reconstruction process, we used this method for resolutions of 2.0, 1.0, 0.8, and 0.5 mas. We compared the resulting reconstructions to the source image, convolved to a resolution corresponding to the reconstruction: we coalligned the reconstruction to match the position of the convolved source image and used the l1-norm as our metric for comparison, which Gomes et al. (2017) found to be the best metric for assessing the quality of a reconstructed image.

Our comparison found that the best reconstructed images came from the reconstructions at 0.5 mas resolution, which is a significantly higher resolution than expected from the maximum projected baseline of the observations. Thus, we opted to use the regularizers and strengths of that reconstruction, but also to reconstruct using the parameters of the next best reconstruction at a lower resolution, in this case at 1.0 mas. For the 0.5 mas resolution, we ran SQUEEZE using the same parameters as the best reconstruction of the simulated data: a 32 × 32 pixel grid, using both total variation and Laplacian regularizers, a mask of a 11 mas diameter uniform disk, an initial image of a 10.5 mas diameter uniform disk, 2000 elements, and 500 iterations. However, in this case we use 25 chains in order to better account for artifacts due to the reconstruction process. Once more we produce a single average and standard deviation image from the result of each chain. Reduced χ² = 1.53 for the November dataset and χ² = 2.08 for the December image are obtained. We follow the same procedure for the 1.0 mas reconstruction, the differences being that our image was 16 × 16 pixels and that we used total variation and the l1-norm of the à trous wavelet transform (Holschneider et al. 1989). We find the reduced χ² = 1.78 for November and χ² = 2.21 for December.

In order to assess the reliability of the reconstruction process, we also used the Multi-aperture Image Reconstruction Algorithm (MIRA) to reconstruct images (Thiébaut 2008). In contrast to SQUEEZE, which uses MCMC minimization, MIRA uses gradient descent. To reconstruct our images, we used a 16 × 16 grid with a pixel scale of 1.0 mas pixel⁻¹ and a 32 × 32 grid with a pixel scale of 0.5 mas pixel⁻¹, the same scales as used with SQUEEZE. We used the same initial images and masks as we did with SQUEEZE. Starting with the initial image, we ran MIRA with a smoothness regularizer for 200 evaluations, then used the resulting centered image as an initial image for a reconstruction using a maximum entropy regularizer with 500 evaluations. We then repeated the process using the resulting centered image as a starting point. In each case, we decreased the regularization strength, μ in each iteration for each regularizer. The resulting images from this method (Fig. C.2) share characteristics, including position and size of bright features, with those obtained using SQUEEZE. From now on, we only consider the 0.5 mas resolution SQUEEZE images for the analysis.

Figure 4 represents the resulting mean images for the 0.5 mas resolution for both epochs. Differences are noticeable between the two epochs. As the (u, v) coverage is marginally different (Fig. 1), and considering the standard deviation of the images (Fig. C.1), we suggest thatthe differences come from a temporal evolution of the stellar surface and are not a residual of the imaging process, as predicted by Chiavassa et al. (2011b). The inhomogeneities are discussed further in Sect. 5, and are compared with the best numerical model reproducing the observations.

Fig. 4 Mean SQUEEZE image reconstruction of CE Tau for a pixel size of 0.5 mas. The top (bottom, respectively) image corresponds to the November (December, respectively) dataset. North is up and east is left. The white ellipse in the bottom left corner represents the main lobe of the synthesized beam.

4.3 Comparison with 3D RHD simulations

Contrast variations of the NIR photosphere of RSG stars have been best explained by the presence of convective cells. Three-dimensional radiative hydrodynamics (RHD) simulations have been used to interpret the interferometric signals of several stars. For Betelgeuse, the comparison was successful in both the H and the K bands (Chiavassa et al. 2010; Montargès et al. 2014). However, more recent observations in the H band have shown that at different epochs on Betelgeuse (Montargès et al. 2016) or at an unprecedented angular resolution on Antares (Montargès et al. 2017), these numerical models were unable to reproduce the departure from classical disk models. In the case of Betelgeuse, it could be the consequence of a change in the convective regime of the star (Kervella et al. 2018) that may have lead to material ejection observed by VLT/SPHERE (Kervella et al. 2016). Because it is lacking from previous observations, it is impossible to conclude on Antares. Alternatively, we note that older observations were in agreement with the simulations: the new instrumentation available in interferometry (better angular resolution, four-telescope u, v coverage) may bias the interpretation.

To assess the compatibility between our PIONIER data of CE Tau and the convective models, we used a numerical simulation produced by the COnservative COde for the COmputation of COmpressible COnvection in a BOx of L Dimensions, L = 2, 3 (CO⁵BOLD, Freytag et al. 2012). We used the simulation st35gm03n13 (Chiavassa et al. 2011a) whose characteristics are presented in Table 4. This model does not include rotation or a magnetic field. Montargès et al. (2017) checked that at the 15th lobe of the visibility function, these simulations are not affected by numerical artifacts. We recall that our VLTI/PIONIER data on CE Tau probe only the first three lobes of the visibility function.

Hundreds of temporal snapshots are computed. Intensity images are computed in the PIONIER spectral channels using the 3D pure local thermodynamic equilibrium (LTE) radiative transfer code Optim3D (Chiavassa et al. 2009). The images are scaled to the angular diameter of CE Tau for each epoch. To account for the unknown orientation on the plane of the sky, each image is rotated around its center. We used 36 angle positions between 0° and 180°. Finally, interferometric observables were derived using a fast Fourier transform (FFT) algorithm.

The observables associated with the grid of temporal snapshots and rotation angles were compared to the squared visibilities of our data. The aim of this procedure is not to find a snapshot reproducing the configuration of the convective pattern of CE Tau exactly, but to find one that reproduces its characteristics (size of the cells, contrast). Matching a configuration pattern would require a much greater number of snapshots, due to the statistics of the convective configurations, and would not provide a better physical interpretation. For both epochs, temporal snapshot #091 was the best match, although the best orientation on the plane of the sky is changing. As shown by Montargès et al. (2017), the angular diameter of the intensity images of the simulations has a strong impact on the quality of the match. Therefore, for the best snapshot, we derive intensity images for each rotation angle at different angular diameters: from 9.59 to 10.59 mas for the November dataset and from 9.68 to 10.68 mas for the December dataset by steps of 0.05 mas. The best matches are obtained for an angular diameter of 10.24 ± 0.05 mas (respectively, 9.98 ± 0.05 mas) and a position angle of 0° for both epochs, and give a ${\tilde{\chi }}^2$ of 98.2 (respectively, 118.6) for the November dataset (respectively, December). Our observations of CE Tau are not well matched by 3D RHD simulations.

5 Discussion

To better compare the SQUEEZE reconstructed images with the 3D RHD simulations, we use the method of Wittkowski et al. (2017) and derive the contrast δI_rms∕⟨I⟩ defined in Tremblay et al. (2013). However, as we have only two epochs of observations and one simulation snapshot, we discard the temporal average and compare the elementary elements. Additionally, we correct for the limb-darkening effect by dividing each image by a best fit LDD intensity model. The resulting contrast is represented in Fig. 6 as a function of the upper radius considered. Near the limb of the star, this contrast is biased by the limb-darkening (LD) effect: small errors in the LD modeling can be interpreted as important fluctuations due to the low intensity at the edge of the stellar disk. Therefore, we will only consider the contrast below a radial cut of 4 mas. We also impose a lower cut of 0.5 mas to exclude the central bright pixel of the simulation images.

With a radial cut of 4 mas, the contrast is 5 ± 1% for the November epoch, 6 ± 1% for the December epoch, 16% for the original simulation image, and 12% for the degraded resolution image (with negligible error bars). We also estimated the contrast over a sample of 30 randomly chosen snapshots of the simulation. We obtained an average value of 23 ± 1% for the original images and 16 ± 1% for the degraded resolution images. It appears that CE Tau presents a lower contrast than the 3D RHD simulation, and that the best matching snapshot corresponds to a lower contrast configuration of the convective pattern. This may come from the broadband imaging we did with SQUEEZE by merging all the spectral channels. However, this flattening of the convective pattern should also be observable on the 3D simulation. Therefore, we suggest that this less contrasted convective pattern is real. This can be a consequence of the younger evolutionary stage of CE Tau compared for example with Antares (Sect. 4.1.2 and Ohnaka et al. 2013) or to its stellar parameters, which are quite different from those involved in the convective simulation. In particular, CE Tau is warmer than the simulation and has a smaller stellar radius: its surface gravity is higher than the model. The combination of these differences can lead to a less prominent convective pattern than the numerical model. Producing these simulations requires a lot of computer resources and it is not yet possible to tailor this simulation grid to the needs of the observations of the various individual RSG stars. It is also possible that CE Tau is now experiencing a quieter episode of convective activity. This was previously observed by Cruzalèbes et al. (2015), who derived a centrosymmetry parameter that was higher than the best 3D RHD simulations at that time. Moreover, a change in the convective activity has already been observed on the prototypical RSG Betelgeuse (Montargès et al. 2016; O’Gorman et al. 2017; Kervella et al. 2018).

Additionally, we note that the contrast of the reconstructed images increases monotonically when we increase the radial cut value. This can be the consequences of inhomogeneities lying at the edge rather than at the center of the disk for the images. However, it can also be a consequence of a residual LD effect. For the simulation, the trend is different with a contrast decreasing between the center and the edge of the star before being subjected to the limb effect. This may indicate that the simulation produces more inhomogeneities at the center of the disk.

Fig. 5 Best snapshot of the RHD simulation. Top: intensity image. Center: intensity image degraded to the resolution of the SQUEEZE reconstructed images. Bottom: reconstructed image from a set of interferometric observables matching the (u, v) coverage and noise level of our PIONIER data.

Fig. 6 Contrast of the convective pattern of the best convective simulation snapshot (orange), the same degraded to the reconstruction resolution (blue), the November (green), and December (red) reconstructed images as a function of the upper radius considered.

6 Conclusion

Our VLTI/PIONIER observations, collected at two different epochs, allowed us to derive refined values of the angular diameter of CE Tau. Using archival photometric measurements we were able to derive updated basic stellar parameters for this star, as well as its mass and its age using evolutionary models. These values allowed us to pick up a 3D RHD simulation and the interferometric observations constrained the temporal snapshots and their rotation angle. We were able to compare the best intensity image of the simulation with reconstructed images for our two epochs of observations. The simulation presents a higher radial contrast variation than the images. We suggest thatit may be a consequence of the lower effective temperature of the simulation compared with CE Tau. As the observations on individual RSG stars are becoming more and more precise and numerous, such models become mandatory in order to properly identify the convective activity of these stars. As previous interferometric observationsshowed a vigorous convective activity, we also suggest that CE Tau is undergoing a quieter convective episode.

As already mentioned for the RSG star Betelgeuse (see Sect. 1), Aurière et al. (2016) have shown that the presence of bright spots at the stellar surface can be traced by spectropolarimetric observations. We have collected contemporaneous (Nov. to Dec. 2016) spectropolarimetric observations of CE Tau during a large program on the TBL/Narval instrument. Although the PIONIER and the Narval instruments do not look at the same height in the atmosphere, and so probably probe a different convective scale, the combination of the two techniques can improve our understanding of the convective surface of the star. Comparison between interferometric and spectropolarimetric observations of CE Tau for the period of November–December 2016 will be presented in a forthcoming paper (Tessore et al., in prep.).

Studying the full temporal evolution of the observed convection would require more than two months of coverage. However, these observations represent the basis of future temporal monitoring of the convective pattern of RSG stars. This information is important for assessing the physical recipe of numerical simulations that aim to understand the atmosphere of RSG stars and to explain the origin of their mass loss.

Movie

Animated gif of the two epochs Access here

Acknowledgements

We would like to thank the anonymous referee whose suggestions and comments led to improvements in this article. We are grateful to ESO’s Director-General, Prof. Tim de Zeeuw, for the allocation of observing time to our program, and to the Paranal Observatory team for the successful execution of the observations. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant agreement No. 665501 with the research Foundation Flanders (FWO) ([PEGASUS]² Marie Curie fellowship 12U2717N awarded to M.M.). A.L. acknowledges financial support from “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, France. F.B. acknowledges funding from NSF awards 1445935 and 1616483. We acknowledge with thanks the variable star observations from the AAVSO International Database that were contributed by observers worldwide and used in this research. We used the SIMBAD and VIZIER databases at the CDS, Strasbourg (France)³, and NASA’s Astrophysics Data System Bibliographic Services. This research has made use of the Jean-Marie Mariotti Center’s Aspro⁴ service, and of the SearchCal service⁵ (co-developed by FIZEAU and LAOG/IPAG). This research made use of Matplotlib (Hunter 2007), Astropy⁶, a community-developed core Python package for Astronomy (Astropy Collaboration et el. 2013), and Uncertainties⁷: a Python package for calculations with uncertainties.

Appendix A Observation log

Appendix B χ² maps for LDD and Gaussian spots models

Fig. B.1 χ2 maps for the single-spot model developed in Sect. 4.1. Left: November dataset. Right: December dataset.

Appendix C Image reconstruction: additional figures

C.1 SQUEEZE reconstruction

In addition to the mean reconstructed images presented in Sect. 4.2, we present in Fig. C.1 the mean images plus and minus the standard deviation.

Fig. C.1 Mean reconstructed image (see Sect. 4.2) minus the standard deviation (left column) and plus the standard deviation (right column). The top row corresponds to the November dataset and the bottom row to the December dataset.

C.2 MIRA reconstruction

Figure C.2 shows the reconstructed image of CE Tau obtained with MIRA. For more details on the reconstruction process with this algorithm, see Sect. 4.2.

Fig. C.2 MIRA reconstructed images of CE Tau at a 0.5 mas resolution. The November image is at the top and the December image at the bottom.

References

All Tables

All Figures

	Fig. 1 (u, v) coverage of our three epochs of VLTI/PIONIER observations of CE Tau. Left: 14 November 2016. Center: 22 November 2016. Right: 23 December 2016. The underlying image corresponds to the best fit LDD model of Sect. 3
In the text

	Fig. 2 Squared visibility in the first two lobes measured by VLTI/PIONIER on CE Tau at 1.62 μm. The purple squares correspond to the November dataset. The orange dots correspond to the December dataset. The continuous curves are the best fit LDD power-law models for both epochs.
In the text

	Fig. 3 Hertzsprung–Russell diagram with the position of CE Tau marked by the red cross. The continuous (respectively, dashed) lines correspond to the non-rotational (respectively, rotational) evolutionary tracks of Ekström et al. (2012) for solar metallicity abundances (Z = 0.014).
In the text

	Fig. 4 Mean SQUEEZE image reconstruction of CE Tau for a pixel size of 0.5 mas. The top (bottom, respectively) image corresponds to the November (December, respectively) dataset. North is up and east is left. The white ellipse in the bottom left corner represents the main lobe of the synthesized beam.
In the text

	Fig. 5 Best snapshot of the RHD simulation. Top: intensity image. Center: intensity image degraded to the resolution of the SQUEEZE reconstructed images. Bottom: reconstructed image from a set of interferometric observables matching the (u, v) coverage and noise level of our PIONIER data.
In the text

	Fig. 6 Contrast of the convective pattern of the best convective simulation snapshot (orange), the same degraded to the reconstruction resolution (blue), the November (green), and December (red) reconstructed images as a function of the upper radius considered.
In the text

	Fig. B.1 χ2 maps for the single-spot model developed in Sect. 4.1. Left: November dataset. Right: December dataset.
In the text

	Fig. C.1 Mean reconstructed image (see Sect. 4.2) minus the standard deviation (left column) and plus the standard deviation (right column). The top row corresponds to the November dataset and the bottom row to the December dataset.
In the text

	Fig. C.2 MIRA reconstructed images of CE Tau at a 0.5 mas resolution. The November image is at the top and the December image at the bottom.
In the text