© ESO 2019

1 Introduction

Red supergiants (RSGs) represent the late stage of the evolution of massive (>8 M_⊙) stars before they explode as type II supernovae. These stars are characterized by high luminosities (~10⁵ L_⊙), low effective temperatures (ranging between 3450 and 4300 K), and low surface gravities (ranging between log g = −0.5 and 0.5); such low gravities result in large radii on the order of 10³ R_⊙ (Levesque et al. 2005; Josselin & Plez 2007).

One of the distinctive properties of RSG stars is their high mass-loss rates (Ṁ ~ 10⁻⁷−10⁻⁴ M_⊙ yr⁻¹; Meynet et al. 2015), which place them among the major contributors to the chemical enrichment of the interstellar medium. However, the driving mechanism of the mass-loss process in these stars is still poorly understood. Josselin & Plez (2007) propose that convective motions, which decrease theeffective gravity, in the atmosphere of RSG stars combined with the radiative pressure on molecular lines could initiate the mass loss. This scenario was later supported by Arroyo-Torres et al. (2015) who found that the extension of the molecular layers increases with both increasing luminosity and with decreasing surface gravity in a sample of RSG stars. An additional mechanism that may contribute to the mass-loss process is the magnetic field, which was recently detected around RSGs (for example, Aurière et al. 2010; Tessore et al. 2017).

Another relevant property of RSG stars is their photometric variability. Kiss et al. (2006) analyzed visual light curves of a sample of RSG stars. These light curves are characterized by V -band amplitudes between 1 and 4 mag. In general, various timescales are involved in their (pseudo-)periodicity. However, the common behavior is represented by two main photometric periods: a short one, on the order of a few hundred days; and a long one, on the order of a few thousand days. The short photometric periods were attributed to atmospheric pulsations in the fundamental or low-overtone modes, while magnetic activity, binarity, or non-radial gravity modes are sometimes invoked to explain the long periods. In addition, according to Schwarzschild (1975), large convective cells in the atmospheres of RSG stars could induce photometricirregularities.

Though more luminous, RSG stars bear some similarities with asymptotic giant branch (AGB) stars of the Mira type. It is shown by Arroyo-Torres et al. (2015) that the (relative) extension of the molecular layers is comparable in RSG stars and in Mira stars. On the other hand, Mira stars are fundamental-mode pulsators (Wood et al. 1999) and are thus characterized by regular photometric variations. Moreover, the V -band amplitudes in Mira stars are between 2.5 and 7.0 mag, which are, in general, larger than those of RSG stars. The mass-loss process is also different in Mira stars, where not only convection, but mainly the combination of pulsations and radiation pressure on dust moves the stellar material outward. Magnetic fields may also contribute to the mass loss (Lèbre et al. 2014).

The present paper aims to relate the atmospheric dynamics in the RSG star μ Cep to its photometric variability with the help of high-resolution observations and state-of-the-art three-dimensional (3D) stellar convection simulations of RSG stars. For this purpose, the recently developed tomographic method, which aims to probe velocity fields at different depths in the stellar atmosphere, is used.

The paper is structured as follows. Section 2 introduces the tomographic technique as well as our method used to infer the effective temperature. The tomographic method was then applied to the RSG star μ Cep (Sect. 3) as well as to snapshots from 3D radiation-hydrodynamics (RHD) simulations (Sect. 4). Section 5 summarizes our conclusions.

2 Methodology

In order to constrain and characterize the atmospheric motions in μ Cep, we aim to determine its spatially-resolved velocity field along with its effective temperature as a function of time. Both are then compared to those provided by 3D RHD CO5BOLD simulations. For this purpose, different techniques are used. They are described in the following sections.

2.1 Tomography

The derivation of the velocity field is performed using the tomographic technique developed by Alvarez et al. (2001a). The tomographic method aims to recover the temporarily-resolved line-of-sight velocity distribution as a function of depth in the stellar atmosphere (within different optical-depth slices). The method was further developed and validated on 3D RHD CO5BOLD simulations of RSG stars by Kravchenko et al. (2018). In short, the tomographic technique is based on sorting spectral lines according to their formation depths as indicated by the maximum of the contribution function (cf., see Kravchenko et al. 2018). The formation depth is expressed on an optical depth scale computed at the reference wavelength λ = 5000 Å. The atmosphere is split into different slices corresponding to specific optical depth ranges. For each slice, a spectral mask is then constructed,which contains the wavelengths of lines forming in the corresponding range of optical depths. For a better wavelength precision, only atomic (and not molecular) lines are kept in masks. The masks are then cross-correlated with either observed or synthetic stellar spectra. The resulting cross-correlation function (CCF) profile reflects the average shape and radial velocity (RV) of lines forming in a given optical depth range.

2.1.1 Construction of masks

A set of five tomographic masks was constructed from a one-dimensional (1D) model atmosphere computed with the MARCS code (Gustafsson et al. 2008). The parameters of the model are T_eff = 3400 K, log g = −0.4, mass = 5 M_⊙, microturbulence velocity 2 km s⁻¹, and solar metallicity. The selected stellar parameters fall in the parameter range of RSGs (in terms of temperature and surface gravity; Levesque et al. 2005) as well as of current 3D simulations (in terms of temperature, surface gravity, andmass; Chiavassa et al. 2011). In the following, we applied this set of five masks to the RSG μ Cep, as well as toa 3D simulation. The properties of the masks, that is, the optical depth boundaries and number of lines per mask, are shown in Table 1.

This “new” set of masks is different from the “old” one used by Kravchenko et al. (2018) in terms of the number of masks and the corresponding optical-depth ranges. The old set of masks was built from a 1D MARCS model atmosphere with the same stellar parameters as the one used in the present work but with a different atmospheric extension. The atmospheric extension may be expressed by the Rosseland optical depth τ_Ross (i.e., the optical depth derived from the Rosseland mean opacity). In general, MARCS model atmospheres are computed between τ_Ross = 10⁻⁵ and 100. In our case, the constructed model extends further to τ_Ross = 10⁻⁶.

2.1.2 CCF analysis

The tomographic masks obtained as described in the previous section are cross-correlated with μ Cep and 3D snapshot spectra in Sects. 3.2 and 4.1, respectively. Cross-correlation is a powerful technique that allows one to sum the information from hundreds of spectral lines (which share properties defined during the masks construction) into an average profile, thus increasing the signal-to-noise ratio (S/N). The quantity $$p\left(\lambda \right)=\genfrac{}{}{0.1ex}{}{1-\text{CCF}\left(\lambda \right)}{{\int }_0^{\infty }[1-\text{CCF}\left(\lambda \right)]\text{d}\lambda }$$(1)

may be interpreted as a line probability distribution (Lion et al. 2013). Thus, the first (raw) moment of p(λ) (hereafter M₁) allows one to characterize the CCF in terms of the mean RV shift.

In the case of an asymmetric (or double-peaked) CCF, one may derive separately the RV for each CCF component. For this purpose, a fit of the CCF with a single, or multiple, Gaussian profile is performed with the help of the detection of extrema (DOE) tool (Merle et al. 2017). In brief, DOE takes the CCF of a given spectrum as input and returns the number of peaks present in the CCF. To this end, DOE computes the first, second, and third derivatives of the CCF by convolving it with the first, second, and third derivative of a narrow Gaussian kernel, respectively. This technique allows one to differentiate the CCFs and smooth them simultaneously, which avoids the numerical noise that may appear when one does discrete differentiation. Then, DOE looks for ascending zeros of the third derivative, that is, inflection points (or local minima of the second derivative), to disentangle the CCF components. The CCF is then fit by a multi-Gaussian function over a small range around the local maxima and/or inflection points in order to obtain the velocity of the identified peaks. Finally, for CCFs with a Gaussian profile, one may assess the full width at half maximum of the CCF as $2.355\sqrt{{\text{M}}_2}$, where M₂ is the second (central) moment of p(λ), that is, its variance.

2.2 Effective temperature determination

The derivation of the effective temperature from observed and synthetic spectra is performed following the method described in Van Eck et al. (2017). This method is based on the computation of band-strength indices $$B=1-\frac{({\lambda }_{C_f}-{\lambda }_{C_i})}{({\lambda }_{B_f}-{\lambda }_{B_i})}\frac{{\displaystyle {\int }_{{\lambda }_{B_i}}^{{\lambda }_{B_f}}{F}_{\lambda }}\text{d}\lambda }{{\displaystyle {\int }_{{\lambda }_{C_i}}^{{\lambda }_{C_f}}{F}_{\lambda }}\text{d}\lambda },$$(2)

where F_λ is the observed or synthetic spectrum, ${\lambda }_{C_i}$ is the lower boundary of the local continuum window, ${\lambda }_{C_f}$ is the upper boundary of the local continuum window, ${\lambda }_{B_i}$ is the lower boundary of the band window, and ${\lambda }_{B_f}$ is the upper boundary of the band window.

A subset of five TiO bands (see Table 2) was selected from the sample of Van Eck et al. (2017). The measured band-strength index is then the mean of the single-band indices from Table 2.

To derive the effective temperature, we proceeded as follows. First, band-strength indices were derived for a grid of 1D MARCS model atmospheres with different effective temperatures in the range of 3400–4400 K, log g = −0.5, and solar metallicity. The relation between the effective temperature of 1D model atmospheres and band-strength indices was thus obtained and is shown in Fig. 1. Then, band-strength indices were derived for either observed or synthetic spectra and the corresponding effective temperatures were deduced from the above calibration (see Sects. 3.2 and 4.1 regarding the application to μ Cep and a 3D simulation, respectively).

3 The red supergiant star μ Cep

Our target star μ Cep was selected from the sample studied by Josselin & Plez (2007) and it is one of the largest and brightest RSGs in our Galaxy. Atmospheric parameters of μ Cep are listed in Table 3; they are taken from Josselin & Plez (2007) and Levesque et al. (2005). The photometric periodicity of μ Cep is analyzed by Kiss et al. (2006). The authors computed the power spectrum of the μ Cep visual light curve and deduced two main photometric periods: a short one of about 860 days and a long one of about 4400 days.

de Wit et al. (2008) derived the mass-loss rate of μ Cep from the presence of an extended circumstellar shell at 25 μm and they conclude that it is on the order of a few 10⁻⁷ M_⊙ yr⁻¹. More recently, De Beck et al. (2010), Mauron & Josselin (2011), Shenoy et al. (2016), and Montargès et al. (2019) derived mass-loss rates of a few 10⁻⁶ M_⊙ yr⁻¹ for μ Cep. This value is much smaller than that of other RSGs, where the mass-loss rate may be orders of magnitude larger (Schuster et al. 2006). VY CMa for instance has a mass-loss rate of 3 × 10⁻⁴ M_⊙ yr⁻¹ (Humphreys et al. 2005).

The atmospheric dynamics of μ Cep was previously studied by Josselin & Plez (2007) by applying the tomographic masks of Alvarez et al. (2001a) to a small sample of ELODIE (Baranne et al. 1996) spectra having a spectral resolution of R ~ 40 000 and covering a time span of about 400 days. Josselin & Plez (2007) detected strong asymmetries in the resulting CCF profiles, which are attributed to complex supersonic velocity fields in the μ Cep atmosphere. They observed as well the time-variable upward and downward motions.

Fig. 1 Band-strength index B as function of effective temperature as computed from synthetic spectra of 1D MARCS model atmospheres.

3.1 Observations

The long-term monitoring of μ Cep was performed with the HERMES spectrograph (Raskin et al. 2011) mounted on the 1.2 m Mercator telescope installed at the Roque de los Muchachos Observatory (La Palma, Spain). The spectral resolution of HERMES is about 86 000 for the high-resolution fiber, and spectra cover the range of 3800–9000 Å.

In total, 95 high-resolution spectra were obtained for μ Cep between April 2011 and January 2018, corresponding to a time span of about 2500 days (i.e., ~ 7 yr). The spectra have a typical S/N of ~100 at 5000 Å. All spectra were reduced using the standard HERMES reduction pipeline, as described by Raskin et al. (2011).

Fig. 2 Excerpt of CCF sequence corresponding to time range between Julian date (JD) 2 456 432 and JD 2 456 607 (compared with light curve of μ Cep on bottom panel of Fig. 3). The last three JD digits are listed on the top of each subpanel. Colors correspond to different masks, mask C1 probes the innermost atmospheric layer while mask C5 probes the outermost layer. Vertical lines at 21.4 km s−1 indicate the mean radial velocity, computed from all masks and epochs (see text).

3.2 Results

The μ Cep HERMES spectra were cross-correlated with the set of tomographic masks constructed in Sect. 2.1.1, and the corresponding RVs were derived as explained in Sect. 2.1.2: a fit of the CCFs with single or multiple Gaussian function(s) was performed in order to extract RVs separately for each CCF component. During the cross-correlation process, the spectra were corrected for the Earth’s motion.

As a next step, temperatures were derived for all μ Cep spectra as described in Sect. 2.2. Our derived temperatures are consistent with those obtained by Levesque et al. (2005) and Josselin & Plez (2007) (see Table 3).

A small set of CCFs is shown in Fig. 2. The CCFs associated with the innermost masks (C1–C2) are shallower than those in the outer layers (C4–C5), since the former are probing weaker lines, which are characterized by higher excitation potentials. The CCF widths increase when going from the inner to the outer atmospheric layers. The average CCF width in mask C1 is about 18 km s⁻¹ while it is about 44 km s⁻¹ in mask C5.

The CCF profiles show asymmetries in all masks. Moreover, at several epochs, the CCFs in masks C1, C2, and C3 that probe the innermost atmospheric layers are characterized by an additional blue (or red) component shifted by 10–30 km s⁻¹ with respect to the main peak. Such asymmetries were not observed by Josselin & Plez (2007), which might be a consequence of the lower resolution of their spectra (twice lower than HERMES), the more restricted spectral range (3900–6800 Å, as compared to 3800–8900 Å for HERMES), and/or the different tomographic masks used.

Line doubling is observed in Mira stars (e.g., Alvarez et al. 2001a,b) and is associated with the passage of a shock wave through the atmosphere (known as the Schwarzschild scenario; Schwarzschild 1954). According to this scenario, when the shock front is located far below the line-forming region, the photosphere contains only falling material so that the spectral lines are red-shifted. Then, as the shock wave approaches the photosphere, a blue component appears and strengthens. It is important to note that this occurrence is associated with the rising matter. When the shock wave passes through the photosphere, the line is split in two components. Finally, when the shock wave has passed through the photosphere, all material is rising, and the line is fully blue-shifted at that point. Alvarez et al. (2000) showed that, in Mira stars, the tomographic method clearly unravels the Schwarzschild scenario as a line-doubling (or shock) front progressing toward external layers as time passes. In other words, in Mira stars, the combined temporal and spatial variations of the CCFs clearly reveal the upward motion of the shock front (see as well Figs. 5 and 6 of Jorissen et al. 2016). In contrast, the temporal and spatial evolution of CCFs in Fig. 2 does not follow the Schwarzschild scenario.

In masks C1 and C2, at least three velocity components are observed as follows. The main CCF component is always located around 10–30 km s⁻¹, whereas additional CCF components appear at various variability phases either at 30–45 km s⁻¹ or at 0–10 km s⁻¹ (see also top panel of Fig. 3). Moreover, the additional CCF component remains weak and never reaches the same contrast as the main CCF peak.

Since the center-of-mass (CoM) velocity of μ Cep is unknown, the distinction between rising and falling material is uncertain. As a proxy of the CoM velocity, we adopt the mean RV of 21.4 km s⁻¹ (an average over all masks and all epochs; it is represented by the vertical lines in Fig. 2 and the horizontal lines in Figs. 3 and 4). To evaluate the mean RV, the (first-moment) M₁ velocity was used (see Sect. 2.1.2).

The intensity of a CCF component is directly related to the size of the corresponding emitting surface on the star, as shown in Fig. 12. For instance, if a large fraction of the stellar surface is covered by rising material, a strong blue-shifted peak is observed in the CCF. Thus, asymmetries in the CCFs hint at the presence of a few granules (or alternatively, at the presence of non-spherically-symmetric shock waves) in the atmosphere of μ Cep which, in any case, behaves differently than the atmospheres of Mira stars.

Figure 3 compares the visual light curve extracted from the American Association of Variable Star Observers (AAVSO) database (Kafka 2018) to the Gaussian-fit or M₁ velocities in the different masks as well as to the TiO-band temperatures. First, Fig. 3 reveals that the TiO-band temperaturevariation is mostly in phase with the light variation. This agreement, in turn, confirms the reliability of our temperature-derivation method (see Sect. 2.2). Second, the maximum peak-to-peak amplitude of the M₁ velocity variations (middle panel of Fig. 3) is similar for all masks and amounts to about 10 km s⁻¹. The corresponding standard deviations lie in the range of 2.0–2.5 km s⁻¹. This situation contrasts with the results of the 3D RHD simulations discussed in Sect. 4.1 where the M₁ velocity amplitudes are highest in the outermost mask. An interesting result is that the RV amplitudes in μ Cep are smaller than those typical of Mira stars by a factor of at least two (Alvarez et al. 2001b). Third, there is a tendency for blue-shifted CCF components to appear in the innermost masks when the stellar brightness increases (top panel of Fig. 3). Fourth, the RVs in all masks vary with the same periodicity as the light, and temperature, but with a phase shift: the maxima in the light curve occur at later epochs than the velocity maxima. The phase shift¹ between the M₁ curve and the light curve was derived by cross-correlating them for two pseudo light cycles in Fig. 3. A pseudo light cycle isdefined as one oscillation of the light curve (starting at an arbitrary epoch). Three such pseudo cycles are identified by the gray shaded areas in Fig. 3. The derived phase shifts for the first and third cycles are similar and amount to about 0.15 for mask C5 and 0.20 for mask C2. Thus, the phase lag decreases toward the outer atmospheric layers.

A similar correlation between RV and light variations exists for Mira stars (e.g., R CMi; see Fig. 5 of Lion et al. 2013) and it is characterized by a phase shift of about 0.15–0.2. Interestingly, in their tests of dynamical DARWIN models for M-type AGB stars, Liljegren et al. (2017) find that a phase shift of 0.2 between the luminosity and radius maxima has to be introduced in order to reproduce the characteristic RV curves of Mira stars.

Such a phase lag was also detected for the RSG star Betelgeuse by Gray (2008). It may be represented as a hysteresis loop between the line-depth ratio (V I 6251.83 Å∕Fe I 6252.57 Å) and the mean core velocity of V I, Fe I, and Ti I (6261.11 Å) spectral lines. The line-depth ratio was considered to be a good temperature indicator and was varying mostly in phase with the light curve. According to Gray (2008), the hysteresis loop illustrates the convective turn-over of the material in the stellar atmosphere: first, the rising hot matter reaches upper atmospheric layers, then temperature drops as the matter moves horizontally and finally matter falls and cools down. However, the above explanation cannot account for observable time-dependent effects if it relates to purely stationary convection. Therefore, in Sects. 4.2–4.4 below, we present a more detailed discussion of the possible origin of the hysteresis loop, searching for mechanisms introducing non-stationary effects in convection.

Figure 4 displays the M₁ velocity as a function of the TiO-band temperature for the three pseudo light cycles of μ Cep defined above and depicted as shaded areas in Fig. 3. The hysteresis loops of μ Cep in Fig. 4 turn counter-clockwise as do those of Betelgeuse (Gray 2008). Moreover, thanks to the tomographic method, the hysteresis loops in μ Cep are now spatially resolved, in the sense that their properties are now probed as a function of depths in the atmosphere.

The detailed inspection of Fig. 4 and Table 4 reveals that the hysteresis loops of μ Cep are characterized by different timescales, as are the corresponding light cycles. Furthermore, there is a tendency for the RV range of the hysteresis loops of the first and third pseudo-cycles to decrease outward, going from mask C2 to mask C5 (mask C1 is rather noisy and does not follow this trend).

The characteristic timescale of the third pseudo light cycle (Table 4) closely matches the short photometric period of μ Cep (860 days; Kiss et al. 2006). The timescales of the hysteresis loops for the first and second cycles are shorter by factors of 1.5 and 2.7, respectively. A similar match between the timescales of the hysteresis loops and the light variations was reported by Gray (2008) for Betelgeuse. The light curve of Betelgeuse is rather regular with a period of about 400 days (Kiss et al. 2006). The radial-velocity and temperature ranges of the hysteresis loops of Betelgeuse and μ Cep are compared inTable 4 and appear to be similar. This similarity suggests that the same physical process is at work in these objects to trigger the velocity and photometric variations. In order to identify this process, we compare in the next section the results obtained for μ Cep with those provided by the 3D RHD CO5BOLD simulation of a RSG star atmosphere.

Fig. 3 Top panel: RVs derived from fit of μ Cep CCFs with Gaussian functions. RVs from the main CCF component are connected. Only CCF components with depths ≤ 2σ, where σ is the standard deviation in the flat part of the CCFs (measuring the “correlation noise”), are kept. Colors correspond to different masks. Middle panel: M1 velocity of CCFs in all masks. The color coding is the same as in the top panel. Bottom panel: AAVSO visual light curve (gray crosses). Black squares correspond to TiO-band temperatures. Vertical lines in all panels indicate times of maximum light and reveal the phase shift between the light (or temperature) and RV variations, the light maximum occurs after the maximum velocity. Horizontal lines in top and middle panels indicate the mean RV. Gray areas define three pseudo light cycles, which correspond to the hysteresis loops in Fig. 4.

Fig. 4 From left to right: hysteresis loops between TiO-band temperature and M1 velocity for three pseudo light cycles of μ Cep (identified as shaded areas in Fig. 3). From top to bottom: M1 velocity measured in different masks, from C1 probing innermost atmospheric layer, to C5 probing outermost layer. Horizontal lines in all panels indicate the mean RV from all masks. The arrow indicates the direction of evolution along the hysteresis loops.

4 3D radiative-hydrodynamics simulations and detailed radiative transfer

Following Kravchenko et al. (2018), we used the same 3D RHD simulation of an RSG star computed with the CO5BOLD code (Freytag et al. 2012), where the combined compressible hydrodynamics and non-local radiative-transfer equations are solved on a Cartesian grid. The code takes into account molecular opacities; however, the radiation transport is treated in a “gray” approximation and ignores radiation pressure and dust opacities. The model geometry is of the kind referred to as “star-in-a-box”. It is assumed that solar abundances are appropriate for RSG stars. The 3D simulations are time-dependent, characterized by realistic input physics, and reproduce the effects of convection and non-radial waves (Chiavassa et al. 2011). Currently, the 3D simulations do not include a radiative-driven wind. However, in RSGs, the wind velocity that can sustain the mass loss must be very low just above the photosphere due to the high densities there (Ohnaka et al. 2017). Therefore, we do not expect that this missing ingredient jeopardizes the comparison with observations performed in the remainder of this paper.

The basic parameters of the 3D simulation are listed in Table 5. Due to the limited number of currently available 3D simulations of RSG stars, the parameters of the selected 3D simulation differ from those of μ Cep, especially in terms of stellar mass. The computation of a 3D simulation for a 25 M_⊙ star would require many more grid points than a simulation of a 5 M_⊙ star, which is already computationally demanding². Nevertheless, this 3D simulation reproduces μ Cep in terms of surface gravity. Moreover, as is shown later in this section, the dynamical picture in a 5 M_⊙ simulation is very different from that of a lower-mass simulation, which represents a pulsating AGB star. The former shows the absence of global shocks and produces photometric and spectroscopic signatures similar to those observed in μ Cep.

The following consequences are expected from the mismatch between the stellar parameters of μ Cep and the 3D simulation. First, there could be an impact on the granule size x_g, which depends on the pressure scale height (see Paladini et al. 2018, and references therein): $$H_{\text{P}}=\frac{\Re T_{\text{eff}}}{\mu g}$$(3)

where g is the surface gravity, ℜ is the gas constant, and μ is the mean molecular weight. The pressure scale height differs by only 2.5% when considering gravity and temperature from either the 3D simulation or from μ Cep (Josselin & Plez 2007). Since the granule size x_g is roughly proportional to H_P (see Paladini et al. 2018, and references therein), the same (negligible) relative difference ensues for the granule size. However, this means that the ratio of granule size to stellar radius is too large in the simulations since the μ Cep radius is much larger than the simulation radius (compare Tables 3 and 5). Second, as the convective efficiency and velocities depend on the effective temperature, convective velocities may be expected to be slightly higher in μ Cep than in the 3D simulation since the latter has a lower effective temperature. Third, opacities may slightly differ between μ Cep and the 3D simulation. Finally, the mass-loss rate may be different in 5 M_⊙ and 25 M_⊙ stars. Fortunately, as we discussed in Sect. 3, the mass-loss rate of μ Cep is not as extreme as that observed in other massive RSGs; so the mass mismatch between the model and actual μ Cep should not be an issue in that respect.

Figure 5 displays the evolution of the bolometric magnitude over a time span of 11 yr covered by the simulations. This synthetic light curve resembles those typical of RSG stars (for example, see Kiss et al. 2006, and Fig. 3) in terms of periodicity, on the order of a few hundreds days, and irregular pattern. Because the observed light curve is given in the V filter (Fig. 3), whereas the 3D RHD simulations generally provide bolometric light curves (Fig. 5), Fig. 6 compares the visual and bolometric light curves for a restricted time range of the 3D RHD simulation chosen arbitrarily³ and shown as a gray-shaded area in Fig. 5, thus allowing one to estimate the bolometric correction. The absolute visual magnitudes were computed by integrating 3D snapshot spectra in the V band using the transmission curve from Bessell (1990). The resulting bolometric corrections are comparable to those of Levesque et al. (2005) determined from the MARCS model atmospheres of RSG stars. Furthermore, the visual light curve in Fig. 6 resembles those typical of RSG stars in terms of amplitude (~1 mag).

Figure 7 shows the pressure scale height H_P, which is spatially averaged over spherical shells as explained in Chiavassa et al. (2009), as a function of geometrical depth. It reveals significantly different values for H_P above and below the stellar radius (corresponding to τ_Ross = 1). This is discussed in Sect. 4.3.

The location of each tomographic mask with respect to the stellar radius is shown in Fig. 7 as colored lines. The position of each mask is defined by the average radial distance from the stellar center for grid points belonging to the considered mask. The averaging process was performed first over a given snapshot and then over time. The vertical bands in Fig. 7 are one standard deviation fluctuations with respect to the time average. Figure 7 shows that, as expected, the tomographic masks probe distinct atmospheric layers. It must be noted that mask C5 is located at just over 900 R_⊙ from the stellar center, whereas the 3D simulation cube has a half-size of only 800 R_⊙. Mask C5 is therefore partially out of the simulation box, and is enclosed only within its corners, which lie 1150 R_⊙ away from the stellar center (see Fig. 8).

We used the pure local thermodynamic equilibrium (LTE) code Optim3D (Chiavassa et al. 2009) to compute synthetic spectra for the time span of the 3D simulation corresponding to the gray-shaded area in Fig. 5. The code computes the radiative transfer in detail using pre-tabulated extinction coefficients as a function of temperature, density, and wavelength for the solar composition (Asplund et al. 2009). They were constructed with no micro-turbulence broadening, and the temperature and density distributions were optimised to cover the values encountered in the outer layers of the RHD simulations. The code takes into account the Doppler shifts caused by the convective motions. The wavelength range of the computed spectra goes from 3700 to 8900 Å, and the spectral resolution is R = 85 000, which mimics the resolution of the HERMES spectrograph (see Sect. 3). A few issues associated with the numerical resolution of the 3D simulation were encountered and resolved as explained in Appendix A.

Fig. 5 Bolometric magnitude as function of time for 3D RHD simulation. The shaded area corresponds to the subset of snapshots analyzed in Sect. 4.1.

Fig. 6 Absolute visual magnitude (red dots and left scale) and bolometric magnitude (black crosses and right scale) variations during restricted time span of 3D RHD simulation (gray-shaded area in Fig. 5). The absolute visual magnitude is computed by integrating 3D snapshot spectra in the V band (using the appropriate transmission curve).

Fig. 7 Pressure scale height (HP) profiles as function of shell radius for set of 3D snapshots. Vertical lines indicate the radial geometrical distances from the stellar center (averaged over a given snapshot and over time) corresponding to the tomographic masks C1–C5, along with the locations corresponding to τRoss = 1 (denoted as R*) and τ0 = 20 000. The bands are standard deviations of the averages over time.

Fig. 8 Line-of-sight velocity structures for cut through center of four different 3D snapshots (for observer located at bottom of each panel). Black shaded areas on velocity maps correspond to regions probed by mask C5.

4.1 Results

The synthetic spectra obtained as described in the previous section were cross-correlated with the tomographic masks obtained in Sect. 2.1.1. An example of resulting CCFs is shown in Fig. 9. The CCFs appear asymmetric in all masks. Moreover, at several epochs, the CCFs are characterized by an additional blue, or red, component shifted by ~ 10 km s⁻¹ with respect to the main peak. A similar behavior was observed in μ Cep CCFs for masks C1–C3 (Fig. 2). But unlike in μ Cep, the simulated CCFs show hints of double-peaked profiles in the outermost masks C4 and C5. As for μ Cep, the simulated CCFs do not follow the Schwarzschild scenario.

Like those of μ Cep, the simulated CCFs associated with the innermost masks are shallower than those associated with the outer layers. The width of CCF profiles increases going from the inner to the outer atmospheric layers. The average CCF widths in masks C1 and C5 are ~10 and ~ 16 km s⁻¹, respectively, which is a factor of 2–2.5 smaller than those of μ Cep. This indicates that the turbulence velocity fields reached by current 3D simulations are still too low with respect to those observed in real RSG stars. This could result from the following shortcomings of current 3D simulations. First, the numerical resolution (i.e., the number of grid points per pressure scale height) of the 3D simulations is still rather low. The size of each grid point is about 4 R_⊙, which introduces difficulties in resolving the convective structures and the shocks that impact them. Second, the effective temperature of the model used in the 3D simulation is lower than that of μ Cep (compare Tables 3 and 5) and, as noted in Sect. 4, the amplitude of the turbulent velocity field is related to the effective temperature. Finally, the atmospheric extension of the current 3D simulations is still too low (Arroyo-Torres et al. 2015). Larger and more diluted photospheres would encompass higher Mach numbers.

For all simulated CCFs, the corresponding radial velocities were derived as described in Sect. 2.1.2, and the TiO-band temperatures were computed as explained in Sect. 2.2. Following Chiavassa et al. (2009), the effective temperature T_eff of a 3D snapshot is defined as the temperature T at radius R for which L(R)∕(4πR²) = σT⁴ (where σ is the Stefan–Boltzmann constant and L is the luminosity). The bottom panel of Fig. 10 shows T_eff and bolometric magnitude variations for the 3D snapshots. The T_eff amplitude ofvariation explains the M_bol amplitude, and they are mostly in phase. The TiO-band temperature variations are twice as large as the T_eff variations (third panel of Fig. 10). The TiO-band temperature follows the variations of T_eff but lags behind. The TiO temperature is based on TiO-band strength and is, thus, very sensitive to the upper photospheric temperature stratification. The TiO temperature appears, however, as a good proxy of T_eff (and M_bol, and the stellar variability), with possibly a slight phase shift. It is also the only temperature we can derive in μ Cep and use to compare our models with observations.

The evolution of the RVs, M₁ velocities, TiO-band temperatures, and visual magnitudes is displayed in Fig. 10. From this comparison, the following conclusions emerge. As in μ Cep, the temperature variation mimics the light variation (except for the time range between days 1750 and 1900). Unlike in μ Cep, the amplitude of the M₁ velocity variations is larger in the outermost masks C3, C4, and C5 (~10 km s⁻¹ peak-to-peak) than in the innermost masks C1 and C2 (~5 km s⁻¹ peak-to-peak). In terms of standard deviations of the M₁ variations, they increase from 1.5 km s⁻¹ in mask C1 to 3.2 km s⁻¹ in mask C5. These values are similar to those observed in μ Cep (2.0–2.5 km s⁻¹). There is no clear correlation between the episodes of line doubling (top panel of Fig. 10) and the light curve. The appearance of a red-shifted secondary peak is more frequent. Finally, the M₁ velocity varies with the same periodicity as the light (or temperature), but with a phase shift of a few tens of days. Similarly to μ Cep, the maximaof the light curve occur at later epochs than the maxima on the M₁ curves. The derived phase shift between the M₁ curve and the light curve (for the time span corresponding to the pseudo light cycle in Fig. 10) is 0.31 for mask C3 and 0.19 for mask C5. Thus, as in μ Cep, the phase lag decreases from mask C3 to mask C5. Due to the irregular behavior of M₁ in masks C1and C2, it was not possible to derive the corresponding phase shifts.

According to the above results, we conclude that there is a strong similarity between μ Cep and the 3D RHD simulation, especially the existence of a phase shift between the light (or temperature) and the velocity variations, which translates into a hysteresis loop (Fig. 11). The hysteresis-like behavior observed for the 3D RHD simulation resembles that of μ Cep (Fig. 4) and Betelgeuse (Gray 2008). Figure 11 reveals that, unlike in μ Cep, the hysteresis RV range in the 3D RHD simulation decreases toward the interior (see as well Table 4). This is a direct consequence of the larger M₁ amplitude in the outermost masks (middle panel of Fig. 10). Since the CoM velocity of the 3D RHD simulation is known (0 km s⁻¹), the distinction between rising and falling material is possible. It appears that the upper part of the loop corresponds to stationary matter in all masks, whereas the lower part of the loop is increasingly redshifted (falling matter) as one considers layers further up. This behavior is easily accounted for with the help of Figs. 12 and 13. On the upper branch of the hysteresis curve, rising and falling material cover similar areas at the surface, giving rise to a zero radial velocity. On the lower branch of the hysteresis curve, the rising material is still present, but it occupies a smaller surface than the falling material.

For the 3D RHD simulation, similar to μ Cep, the RVand temperature ranges as well as the characteristic timescale of the hysteresis loop were estimated. They are compared to those of μ Cep and Betelgeuse in Table 4 and all share the same order of magnitude⁴. Thus, understanding the origin of the hysteresis loops from the 3D RHD simulation is crucial for correctly interpreting those observed in actual RSG stars. This is the subject of the next section.

Fig. 9 Excerpt of CCF sequence corresponding to time range between days 1990 and 2200 of 3D RHD simulation (see corresponding light curve in bottom panel of Fig. 10), for same masks as in Fig. 2. The day number corresponding to each snapshot is indicated on top of each subpanel. Vertical lines in all panels indicate the 0 km s−1 velocity (i.e., the CoM velocity of the 3D simulation).

Fig. 10 Same as Fig. 3 for snapshots of 3D RHD simulation, except for bottom panel. Bottom panel: Teff (squares) and bolometric magnitude (crosses) variations for 3D snapshots.

Fig. 11 Hysteresis loop between TiO-band temperature and M1 velocity for snapshots from 3D simulation. Colors correspond to different masks, with mask C1 probing the innermost atmospheric layer, and mask C5 the outermost. Horizontal lines in all panels indicate the CoM velocity of the 3D simulation, that is, 0 km s−1. The arrow indicates the direction of the evolution along the hysteresis loop.

4.2 Understanding hysteresis loops from the 3D RHD simulation

The 3D simulations provide an opportunity to follow the evolution of temperature and velocity along the hysteresis loop, and hence to understand its physical origin. Figure 12 displays velocity and temperature maps for different snapshots along the hysteresis loop of mask C4 (Fig. 11). The maps are weighted by the contribution function of a line contributing to the mask C4 (see Kravchenko et al. 2018, for details). The contributionfunction reveals at which depth, and hence velocity and temperature, in the atmosphere spectral lines form.

The velocity maps in Fig. 12 reveal upward and downward motions of matter extending over large portions of the stellar surface. The relative fraction of upward and downward motions is what distinguishes the upper from the lower part of the hysteresis loop; its top part (zero velocity) is characterized by equal surfaces of rising and falling material. The bottom partof the hysteresis loop occurs, as expected, when the stellar surface is covered mostly by downfalling material.

On the other hand, the weighted temperature maps in Fig. 12 can be considered as a proxy to the surface intensity. The maps clearly demonstrate the brightening of those regions of the stellar surface where matter is about to rise. These maps make a clear distinction between the left and right parts of the hysteresis loop; the right part is characterized by the presence of high temperatures at several locations on the stellar surface.

The general trend observed in Fig. 12, as described above, confirms that hysteresis loops reflect the turn-over of material in the stellar atmosphere. The appearance of bright and warm regions at the stellar surface is followed by the rising of material (as seen on the weighted velocity maps) at those same locations, thus accounting for the phase shift observed in Fig. 10.

Figure B.1 displays the weighted temperature and velocity maps along the hysteresis loop of mask C4, which was constructed for a different time range of the 3D simulation (i.e., between days 1500 and 1950). The temperature and velocity maps show the same trend as the one observed in Fig. 12, thus, confirming that the same physical mechanism takes place in any subset of 3D snapshots.

The appearance and disappearance of warm regions on the weighted temperature maps in Fig. 12 are thus responsible for the surface brightness variations. It is important to stress, however, that a crucial property of convection in both AGB and RSG stars, as revealed by the 3D simulations (e.g., Freytag et al. 2017, for AGB stars), is namely the fact that the continuum forms above the top of the convection zone. In both RSG and AGB stars, the deep large-scale convective cells (as displayed on Figs. 13 and 14) are not directly observable. In RSG stars, we see almost down to the top of the convection zone. The continuum-forming layers sit so close to the top of the convection zone that both move together⁵. However, the structures we see are non-stationary surface granules that are affected by acoustic waves or pulsations.

In contrast, in some of the cooler AGB models, we often see much higher layers that are far above the convection zone. Here, the structures can still be shaped by convection because waves have traveled through the top of the convection zone and thus have been shaped by variations in sound speed, density, and velocity. These stochastic shocks generated and shaped by convection transfer the heat through the atmosphere and, in turn, cause the surface brightness variations.

4.3 Surface versus deep convection

According to Fig. 7, the velocity and temperature maps of Fig. 12 describe the convective pattern located well above the bottom of the photosphere since the atmospheric layers probed by the tomographic masks are far up the Rosseland radius, thus marking the stellar surface. The corresponding pressure scale-height values are relatively small (reflecting small-scale surface granulation) with respect to those in the deep convective zone (defined as the region below the Rosseland radius) where large convective cells are located (see Figs. 13 and 14). An increase of the size of convective structures with increasing atmospheric depth is observed as well in convection simulations of main-sequence F-type stars (Kitiashvili et al. 2016).

Figure 14 shows the velocity and temperature maps corresponding to the deep convection zone at τ₀ = 20 000 (see Fig. 7) for the same 3D snapshots as in Fig. 12. The large convective cells are clearly visible as bright granules on the temperature maps. Following the velocity maps, the matter rises through the granules and falls in the intergranular lanes, thus it is in accordance with the classical convection scenario.

The rising material observed in Fig. 12 originates in the deep convective zone and an atmospheric shock develops higher up along the following steps. Freytag et al. (2017) and Liljegren et al. (2018) illustrate this process in similar 3D RHD simulations as follows. Non-stationary convection (e.g., merging downdrafts or other localized events) produces a sound wave in the stellar interior. The sound wave, then, travels through the star until it hits the surface. At the surface, the wave is slowed down and compressed due to the drop in temperature and sound speed. In addition, the amplitude of the sound wave rises due to the decrease in density and, therefore, turns into a shock. A shock, then, propagates all the way from the stellar surface to the outer atmospheric layers (see Fig. 13).

While this happens for the Sun far out in the chromosphere, it happens for RSGs already in the photosphere, since all velocities and also Mach numbers are larger in RSGs (Freytag & Chiavassa 2013). Therefore, the entire RSG spectrum and also the surface brightness are affected, and not only some exotic chromospheric emission lines as in the Sun.

The link between the maps in Figs. 12 and 14 is best illustrated by Fig. 13 (the online version of the paper displays an animated version of this figure), which presents radial-velocity maps in an equatorial plane for the snapshots along the hysteresis loop of mask C4 in Fig. 11. It reveals that convection is structured through the motions going all the way from the center of the star to the surface, which are separated by the downdrafts appearing as fingers on Fig. 13. That figure clearly reveals the structural change occurring between deep and surface convection.

Fig. 12 Top panel: line-of-sight velocity maps (for observer located in front of figure) weighted by CF for snapshots along hysteresis loop of mask C4. Red color corresponds to falling material, and blue color to rising material. Bottom panel: same as top panel for temperature. The arrow indicates the direction of the evolution along the hysteresis loop. The horizontal ticks on the top velocity map show the location of the equatorial plane of 3D snapshots displayed in Fig. 13.

Fig. 13 Velocity maps computed along rays from center of star in an equatorial plane for snapshots along hysteresis loop of mask C4 (online version of paper displays animated version of this figure). This figure displays the large-scale structure of the convection, from the stellar center to the stellar surface. Moreover, it illustrates the appearance of shocks at the stellar surface and their propagation through the atmosphere.

4.4 Convective and acoustic timescales

As an important preamble to this section, one should realize that stationary convection cannot generate time-dependent effects like a hysteresis loop. Its origin must be looked for in time-dependent effects. These could originate from acoustic waves. These waves could disturb the top of the convective pattern, just below the continuum-forming region, and move to the line-forming regions where they modify the flow pattern, as discussed in Sect. 4.3. In this section we investigate whether convective or acoustic timescales may account for the characteristic timescale of the hysteresis loops in RSGs (Table 4), which is of the same order as the observed photometric variations. As explained above (Sect. 4.3), the surface features are probably not genuine convection, but are rather linked to shock waves originating in acoustic waves generated by the deep convection. Therefore, we shall consider both convective and acoustic (i.e., sound-crossing) timescales.

The sound velocity can be expressed as: $$v_{\text{s}}=\sqrt{{\Gamma }_1\frac{\Re T_r}{\mu }},$$(4)

where Γ₁ is the first adiabatic exponent (Γ₁ = 5∕3 for a mono-atomic perfect gas), T_r is the stellar temperature at radial distance r from the stellar center, and μ = 1.3 g mol⁻¹ for the atmosphere of a late-type star (Chiavassa et al. 2011). A single convective timescale is much harder to define, as there is a whole range of them, related to the different convective spatial scales, all the way from the big convective cells to turbulent eddies. Nevertheless, to fix the ideas, one may use the mean lifetime of a turbulence element in the convective zone as expressed by the mixing length theory (MLT) of convection (Böhm-Vitense 1958): $$t_{\text{c}}=\frac{\Lambda }{v_{\text{c}}},$$(5)

where Λ ≃ H_P is the mixing length (i.e., the vertical length traveled by a convective globule before dissolving in its environment) and v_c is the velocity of the convective globule. Following Lamers & Levesque (2017), the convective velocity can be derived from the relation: $$v_{\text{c}}^3=\frac{L_r}{6\pi r^2}{(\frac{{\rho }_r\Re }{\mu })}^{-1}\frac{g_r\Lambda }{T_r},$$(6)

where T_r, L_r, ρ_r, and g_r are the temperature, stellar luminosity, density, and gravity at distance r from the stellar center.

The sound and convective velocities are first derived for the deep interior (r = R∕2, where R is the stellar radius) and then for the surface (r = R), adopting the relevant physical quantities from a snapshot of the 3D simulation. The temperature T_r, density ρ_r, and pressure scale height H_P(r) are computed at r = R∕2 and R by averaging the 3D snapshot variables over spherical shells (as explained in Chiavassa et al. 2009). In the deep interior, Eqs. (4) and (6) predict v_s = 33.8 km s⁻¹ and v_c = 15.2 km s⁻¹, whereas in the outer layers v_s = 6.3 km s⁻¹.

The convective velocity derived from Eq. (6) above may be compared to the value directly extracted from the 3D snapshot. For each grid point of the considered 3D snapshot, we computed the convective velocity as $v_{\text{c}}^2=(v_x^2+v_y^2+v_z^2)$, where v_x, v_y, and v_z are velocities along the x, y, and z directions of the simulation box. The v_c velocity averaged in the range of 0.45–0.55 R amounts to 12.3 km s⁻¹ and agrees well with the value derived from Eq. (6).

In order to derive a convective timescale amidst an extended spectrum of timescales, according to the above caveat, in the deep interior of the 3D simulation, we use Eq. (5) with Λ =H_P(r = R∕2) ≃ 60 R_⊙ (see Fig. 7). The corresponding timescale is on the order of 30 d, that is, it is much too short to account for those involved in the hysteresis loops (Table 4). Further up in the atmosphere, the timescale would become even shorter since H_P decreases outward. Further down in the star where large convective cells are located, the respective timescale would reachan order of years. Thus, the convective timescales cover a wide range from weeks to years, which would comprise the 324-day timescale of the simulated histeresis loop at some specific depth in the star. However, convection is not able to produce a particular, clearly observable feature with that precise timescale.

Instead,the acoustic timescale in the outer layers is a much better choice, assuming the vertical length traveled by a shock or acoustic wave is equal to the distance between the stellar surface (at r = R = 582R_⊙ in Fig. 7) and to the atmospheric layer probed by mask C4 (at r ≃ 800R_⊙ in the same figure). By adopting v_s = 6.3 km s⁻¹ in the outer layers as obtained above, the sound-crossing timescale amounts to 279 d, which is close to the 324 d characteristic time of the hysteresis loop reported in Table 4. This confirms our above mentioned guess that hysteresis loops must have an acoustic origin, although they are not fully independent of convection since the acoustic wave modulates the surface flow that keeps the memory of the underlying convective structure (Freytag et al. 2017).

Fig. 14 Same as Fig. 12, but for τ0 = 20 000. Variations reflected by the hysteresis loop are not observed at this optical depth.

5 Conclusions and future prospects

The present paper applies the tomographic method to the RSG star μ Cep in order to recover its line-of-sight velocity distribution, which is over the stellar disk and within different optical-depth slices, and to relate it to the photometric variations. The observed velocity variations follow the photometric and temperature variations with a phase lag. This phase lag results in hysteresis loops in the temperature–velocity plane that are characterized by timescales of a few hundred days, which are similar to the photometric ones. The same behavior was observed by Gray (2008) for the RSG star Betelgeuse.

The hysteresis loop was also detected in the 3D RHD CO5BOLD simulation of an RSG star atmosphere. The qualitative agreement between the amplitudes of RV and temperature variations as well as of the timescales between observed and simulated hysteresis loops indicates that physical processes related to convection are responsible for the few-hundred-day photometric variations in μ Cep and Betelgeuse.

The timescales of observed and simulated hysteresis loops were compared to theoretical predictions of convective and sound-crossing timescales. It is found that the sound-crossing timescale in the outer layers is of the same order as the hysteresis-loop timescale. This suggests that hysteresis loops are linked to acoustic waves originating from a disturbance in the convective flow below the photosphere and propagating upward to the surface layers where they modulate the convective energy flux.

Perspectives of this work include the application of the tomographic method to an extended sample of RSG stars in order to investigate whether the presence of hysteresis loops is a common feature among them. This study is deferred to a forthcoming paper.

Movie

Movie of Fig. 13 Access here

Acknowledgements

K.K. acknowledges the support of a FRIA (FNRS) fellowship. S.V.E. thanks to Fondation ULB for its support. This work is based on observations obtained with the HERMES spectrograph, which is supported by the Fund for Scientific Research of Flanders (FWO), Belgium, the Research Council of K.U. Leuven, Belgium, the Fonds de la Recherche Scientifique (F.R.S.-FNRS), Belgium, the Royal Observatory of Belgium, the Observatoire de Genève, Switzerland and the Thüringer Landessternwarte Tautenburg, Germany. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research.

Appendix A Numerical resolution of 3D simulations

As was already mentioned in Chiavassa et al. (2009) and Kravchenko et al. (2018), some rays of the 3D simulation are characterized by large differences in the optical depth values between two adjacent grid points near the continuum-forming layers. This, in turn, causes uncertainties in the surface intensity computation. In the case of strong spectral lines, the surface intensity reaches values higher than that of the local continuum and produces spurious emission-like features in line profiles.For example, Fig. A.1 shows continuum-normalized synthetic spectrum around the λ₀ = 4359.64 Å line contributing to mask C5, which probes the outermost atmospheric layer, for a snapshot from the 3D simulation. It is characterized by an “emission” feature with respect to the local continuum in the blue-shifted wing of the spectral line.

Fig. A.1 Continuum-normalized synthetic spectrum around λ0  =  4359.64 Å (around line contributing to mask C5) for 3D snapshot. The wavelength axis was transformed into velocity. The solid line profile was produced using all the rays of the 3D snapshot while the dashed line shows the profile obtained by selecting only rays characterized by a good optical-depth sampling (see main text).

In the meantime, we tried to resolve the issue associated with a poor optical-depth sampling by linearly interpolating the temperature, and the density, on a finer geometrical- or optical-depth grid along 3D snapshot rays. The first approach implied an oversampling of the snapshot box twice (resulting into a 801³ data cube); the corresponding variables (temperature, density, velocity) were interpolated at intermediate grid points. However, this approach allowed us to obtain only two times more data points in the required optical-depth region, but it still lacks enough data for a reliable intensity computation (i.e., with no emission features). The second approach involved an oversampling of the optical depth region between log τ = 0 and 2 by at least factor of five. However, since the true temperature, or density, profile in these regions is unknown, the linear interpolation could not produce a reliable intensity.

An optimal solution was nevertheless found, and it is based on only computing the intensity for rays that are characterized by a good optical-depth sampling near the continuum-forming layer. Several tests were performed in order to find an optimal number of points in the optical-depth range between 1 and 100, which are requiredto obtain the reliable integrated intensity. Our solution requires at least five points between optical-depth values of 1 and 100. It was applied to all 3D snapshots and allowed us to keep at least 90 percent⁶ of rays per snapshot, thus not introducing a significant effect on the overall spectrum (see Fig. A.2). Moreover, the resulting line profile (displayed as the dashed line in Fig. A.1) does not show any artificial emission features anymore. The procedure described above was applied to all of the 3D snapshots analyzed in Sect. 4.

Fig. A.2 Synthetic spectra computed for 3D snapshot with (red) and without (black) removing rays with poor optical-depth sampling. The inset shows a zoom between 5885 and 5900 Å. As a result, most of the emission features are suppressed without affecting the overall spectrum.

Appendix B Additional figure

Fig. B.1 Same as Fig. 12, but for time range between days 1500 and 1950 (Fig. 10).

References

All Tables

All Figures

	Fig. 1 Band-strength index B as function of effective temperature as computed from synthetic spectra of 1D MARCS model atmospheres.
In the text

Fig. 2 Excerpt of CCF sequence corresponding to time range between Julian date (JD) 2 456 432 and JD 2 456 607 (compared with light curve of μ Cep on bottom panel of Fig. 3). The last three JD digits are listed on the top of each subpanel. Colors correspond to different masks, mask C1 probes the innermost atmospheric layer while mask C5 probes the outermost layer. Vertical lines at 21.4 km s−1 indicate the mean radial velocity, computed from all masks and epochs (see text).

In the text

Fig. 3 Top panel: RVs derived from fit of μ Cep CCFs with Gaussian functions. RVs from the main CCF component are connected. Only CCF components with depths ≤ 2σ, where σ is the standard deviation in the flat part of the CCFs (measuring the “correlation noise”), are kept. Colors correspond to different masks. Middle panel: M1 velocity of CCFs in all masks. The color coding is the same as in the top panel. Bottom panel: AAVSO visual light curve (gray crosses). Black squares correspond to TiO-band temperatures. Vertical lines in all panels indicate times of maximum light and reveal the phase shift between the light (or temperature) and RV variations, the light maximum occurs after the maximum velocity. Horizontal lines in top and middle panels indicate the mean RV. Gray areas define three pseudo light cycles, which correspond to the hysteresis loops in Fig. 4.

In the text

Fig. 4 From left to right: hysteresis loops between TiO-band temperature and M1 velocity for three pseudo light cycles of μ Cep (identified as shaded areas in Fig. 3). From top to bottom: M1 velocity measured in different masks, from C1 probing innermost atmospheric layer, to C5 probing outermost layer. Horizontal lines in all panels indicate the mean RV from all masks. The arrow indicates the direction of evolution along the hysteresis loops.

In the text

	Fig. 5 Bolometric magnitude as function of time for 3D RHD simulation. The shaded area corresponds to the subset of snapshots analyzed in Sect. 4.1.
In the text

	Fig. 6 Absolute visual magnitude (red dots and left scale) and bolometric magnitude (black crosses and right scale) variations during restricted time span of 3D RHD simulation (gray-shaded area in Fig. 5). The absolute visual magnitude is computed by integrating 3D snapshot spectra in the V band (using the appropriate transmission curve).
In the text

Fig. 7 Pressure scale height (HP) profiles as function of shell radius for set of 3D snapshots. Vertical lines indicate the radial geometrical distances from the stellar center (averaged over a given snapshot and over time) corresponding to the tomographic masks C1–C5, along with the locations corresponding to τRoss = 1 (denoted as R*) and τ0 = 20 000. The bands are standard deviations of the averages over time.

In the text

	Fig. 8 Line-of-sight velocity structures for cut through center of four different 3D snapshots (for observer located at bottom of each panel). Black shaded areas on velocity maps correspond to regions probed by mask C5.
In the text

	Fig. 9 Excerpt of CCF sequence corresponding to time range between days 1990 and 2200 of 3D RHD simulation (see corresponding light curve in bottom panel of Fig. 10), for same masks as in Fig. 2. The day number corresponding to each snapshot is indicated on top of each subpanel. Vertical lines in all panels indicate the 0 km s−1 velocity (i.e., the CoM velocity of the 3D simulation).
In the text

	Fig. 10 Same as Fig. 3 for snapshots of 3D RHD simulation, except for bottom panel. Bottom panel: Teff (squares) and bolometric magnitude (crosses) variations for 3D snapshots.
In the text

	Fig. 11 Hysteresis loop between TiO-band temperature and M1 velocity for snapshots from 3D simulation. Colors correspond to different masks, with mask C1 probing the innermost atmospheric layer, and mask C5 the outermost. Horizontal lines in all panels indicate the CoM velocity of the 3D simulation, that is, 0 km s−1. The arrow indicates the direction of the evolution along the hysteresis loop.
In the text

Fig. 12 Top panel: line-of-sight velocity maps (for observer located in front of figure) weighted by CF for snapshots along hysteresis loop of mask C4. Red color corresponds to falling material, and blue color to rising material. Bottom panel: same as top panel for temperature. The arrow indicates the direction of the evolution along the hysteresis loop. The horizontal ticks on the top velocity map show the location of the equatorial plane of 3D snapshots displayed in Fig. 13.

In the text

Fig. 13 Velocity maps computed along rays from center of star in an equatorial plane for snapshots along hysteresis loop of mask C4 (online version of paper displays animated version of this figure). This figure displays the large-scale structure of the convection, from the stellar center to the stellar surface. Moreover, it illustrates the appearance of shocks at the stellar surface and their propagation through the atmosphere.

In the text

	Fig. 14 Same as Fig. 12, but for τ0 = 20 000. Variations reflected by the hysteresis loop are not observed at this optical depth.
In the text

	Fig. A.1 Continuum-normalized synthetic spectrum around λ0  =  4359.64 Å (around line contributing to mask C5) for 3D snapshot. The wavelength axis was transformed into velocity. The solid line profile was produced using all the rays of the 3D snapshot while the dashed line shows the profile obtained by selecting only rays characterized by a good optical-depth sampling (see main text).
In the text

	Fig. A.2 Synthetic spectra computed for 3D snapshot with (red) and without (black) removing rays with poor optical-depth sampling. The inset shows a zoom between 5885 and 5900 Å. As a result, most of the emission features are suppressed without affecting the overall spectrum.
In the text

	Fig. B.1 Same as Fig. 12, but for time range between days 1500 and 1950 (Fig. 10).
In the text