Three-dimensional hydrodynamical CO5BOLD model atmospheres of red giant stars VI. First chromosphere model of a late-type giant

TL;DR

This paper presents the first 3D hydrodynamical model of a red giant star's chromosphere, revealing its dynamic shock-driven structure and impact on observable properties, highlighting the need for advanced 3D modeling.

Contribution

It introduces the first 3D hydrodynamical chromosphere model for a late-type giant star, demonstrating its complex dynamics and influence on observable features.

Findings

Chromosphere develops rapidly with shock waves

Peak shock temperatures reach up to 5000 K

3D models better reproduce observable properties

Abstract

Although observational data unequivocally point out to the presence of chromospheres in red giant stars, no attempts have been made so far to model them using 3D hydrodynamical model atmospheres. We therefore compute an exploratory 3D hydrodynamical model atmosphere for a cool red giant in order to study the dynamical and thermodynamic properties of its chromosphere, as well as the influence of the chromosphere on its observable properties. 3D radiation hydrodynamics simulations are carried out with the CO5BOLD model atmosphere code for a star with the atmospheric parameters (Teff=4010 K, log g=1.5, [M/H]=0.0), which are similar to those of the K-type giant star Aldebaran (alpha Tau). ... we compute the emergent continuum intensity maps at different wavelengths, spectral line profiles of Ca II K, the Ca II infrared triplet line at 854.2nm, and H alpha, as well as the spectral energy distribution (SED) of the emergent radiative flux. The initial model quickly develops a dynamical chromosphere characterised by propagating and interacting shock waves. The peak temperatures in the chromospheric shock fronts reach values on the order of up to 5000 K although the shock fronts remain quite narrow. Like for the Sun, the gas temperature distribution in the upper layers is composed of a cool component due to adiabatic cooling in the expanding post-shock regions and a hot component due to shock waves. For this red giant model, the hot component is a rather flat high-temperature tail, which nevertheless affects the resulting average temperatures significantly. The simulations show that the atmospheres of red giant stars are dynamic and intermittent. Consequently, many observable properties cannot be reproduced with one-dimensional static models but demand for advanced 3D HD modelling. Furthermore, including a chromosphere in the models might produce significant contributions to the emergent UV flux.

Equations

$$T_{\mathrm{em}}(z)={⟨{\left(\frac{⟨\kappa \rho T^4{⟩}_{x,y}}{⟨\kappa \rho {⟩}_{x,y}}\right)}^{1/4}⟩}_t$$

Full text

*Methods. *Three-dimensional radiation hydrodynamics simulations are carried out with the CO5BOLD model atmosphere code for a star with the atmospheric parameters ($T_{\mathrm{eff}}\approx 4010$ K, $\log g=1.5$, $[\mathrm{M/H}]=0.0$), which are similar to those of the K-type giant star Aldebaran ($\alpha$ Tau). The computational domain extends from the upper convection zone into the chromosphere ($7.4\geq\log\tau_{\mathrm{Ross}}\geq-12.8$) and covers several granules in each horizontal direction. Using this model atmosphere, we compute the emergent continuum intensity maps at different wavelengths, spectral line profiles of Ca ii K, the Ca ii infrared triplet line at 854.2 nm, and H$\alpha$, as well as the spectral energy distribution (SED) of the emergent radiative flux.

*Results. *The initial model quickly develops a dynamical chromosphere that is characterised by propagating and interacting shock waves. The peak temperatures in the chromospheric shock fronts reach values on the order of up to 5 000 K although the shock fronts remain quite narrow. Like for the Sun, the gas temperature distribution in the upper layers is composed of a cool component due to adiabatic cooling in the expanding post-shock regions and a hot component due to shock waves. For this red giant model, the hot component is a rather flat high-temperature tail, which nevertheless affects the resulting average temperatures significantly.

*Conclusions. *The simulations show that the atmospheres of red giant stars are dynamic and intermittent. Consequently, many observable properties cannot be reproduced with one-dimensional static models but demand for advanced 3D hydrodynamical modelling. Furthermore, including a chromosphere in the models might produce significant contributions to the emergent UV flux.

††offprints: sven.wedemeyer-bohm@astro.uio.no

1 Introduction

The existence of chromospheres around red giant stars is suggested by a variety of observational facts. Mostly, the evidence comes from indirect indicators related to stellar activity, such as, for example, measurements of (variable) flux in the emission cores of Ca II H & K and Mg II h & k lines, asymmetries and shifts of the H$\alpha$ line core (e.g., Meszaros et al. 2008; Vieytes et al. 2011). Observations of red giants in the (sub)millimetre and centimetre ranges also point out to the presence of chromospheres (e.g., Dehaes et al. 2011; Harper et al. 2013). From the theoretical side, however, the properties of red giant chromospheres are still relatively poorly understood. For example, 1D semi-empirical models have been constructed to account for various observable properties of these stars, such as the strengths of chromospheric spectral lines, fluxes in the milimeter wavelength range, and so on (e.g., McMurry 1999). However, these models suffer from various shortcomings. For instance, earlier models were unable to account simultaneously for the required ultraviolet (UV) pumping of carbon monoxide (CO) resonance lines and cool excitation temperatures of these lines in the chromosphere of Aldebaran (McMurry & Jordan 2000). The inability of models to reproduce cool features (e.g., CO lines) and UV lines, led authors to suggest model atmospheres with multiple components and the presence of shock waves in $\alpha$ Tau (McMurry et al. 1999). Detailed multi-component models with a realistic description of time-dependent shock waves require the utilisation of 3D hydrodynamical model atmosphere codes. Unfortunately, to the best of our knowledge, red giant models with chromospheres computed with such codes have not been published so far. Three-dimensional hydrodynamical models are well-suited for studying non-stationary stellar chromospheres and thus may help to better understand the complex interplay between non-stationary physical phenomena (such as, e.g., shock wave activity) that shape their structures.

In this work, we present the results of our first exploratory simulation of the red giant chromosphere performed with the 3D hydrodynamical CO5BOLD model atmosphere package. The atmospheric parameters of the model are similar to those of the K-type giant star Aldebaran, $\alpha$ Tau (see Sect. 2.2 below). Our main goal is to study the dynamical and thermodynamic structure of the chromosphere, and to investigate its influence on the observable properties of this red giant. It must be emphasised that the simulations presented here are only a first step. The included physical effects restrict to what extent the model results can be interpreted. The exploratory model reproduces conditions as they should only be expected in rough approximation for very quiet regions in the low chromosphere. While this is already of interest, $\alpha$ Tau exhibits both C IV and O VI features (Dupree et al. 2005; Robinson et al. 1998). These observations imply significantly higher temperatures than those achieved with shock heating in our 3D hydrodynamical model atmosphere and hints at yet-to-be-included physical (heating) mechanisms.

The paper is structured as follows: after a brief description of the numerical codes and models in Sect. 2, we present the obtained results in Sect. 3, followed by a discussion and conclusions in Sects. 4 and 5.

2 Methodology

2.1 Three-dimensional radiation hydrodynamic simulations

The numerical simulations are carried out with the radiation hydrodynamics code CO5BOLD (Freytag et al. 2012). The hydrodynamic equations and the frequency-dependent radiative transfer equations are solved time-dependently for a fully compressible plasma in a constant gravitational field. The lateral boundaries are periodic. The bottom boundary is open so that material can flow into and out of the computational box. The entropy of the inflowing material is prescribed and sets the effective temperature of the simulated star. The top boundary is transmitting, i.e., material can leave and enter the model at the top. The equation of state is used in form of a look-up table, which is computed in advance under the assumption of thermodynamic equilibrium. For this table, partial ionisation of $\@element{\mathrm{H}}$ and $\@element{\mathrm{He}}$, the formation and dissociation of ${\rm H}_{2}$, and a representative metal are taken into account. The chemical element abundances are adopted from the CIFIST project (Ludwig et al. 2009; Caffau et al. 2011). Non-equilibrium effects like the time-dependence of the hydrogen ionization degree in chromospheres are not considered yet for this first exploratory model due to the high computational costs but should be included in more advanced future detailed models. Therefore, it has yet to be seen how much the atmospheric state for the modelled stellar type would be affected when these effects are taken into account. A long characteristics scheme with multi-group opacities is used for the solution of the radiative transfer equation (see, e.g., Ludwig et al. 1994, for more details on the opacity binning). The opacity data is provided in form of a look-up table with five bins, which are constructed in advance using opacities from the MARCS model atmosphere package (Gustafsson et al. 2008).

2.2 Numerical model

The numerical model analysed here was produced with the radiative hydrodynamics code CO5BOLD (Sect. 2.1 above) in several steps. The initial model atmosphere is based on a snapshot of a well relaxed simulation that was computed using the following atmospheric parameters: $T_{\mathrm{eff}}=4\,000$ K, $\log g=1.5$, and $[\mathrm{M/H}]=0.0$. This snapshot was taken from the CIFIST 3D hydrodynamical model atmosphere grid (Ludwig et al. 2009). The initial model includes the top of the convection zone and the photosphere with a total extent of $4.6\cdot 10^{9}\,\mathrm{m}\times 4.6\cdot 10^{9}\,\mathrm{m}\times 2.8\cdot 10^{9}\,\mathrm{m}\ (x,y,z)$. The size of the initial model was increased to cover more granules by doubling the computational box in each horizontal direction, thus producing $2\,\times\,2$ initially identical quadrants. A multi-scale random velocity field was added to the new model in order to accelerate the diversion of the four quadrants and thus remove the initial symmetry as fast as possible. The sequence was run for one million seconds of stellar simulation time until the initial symmetry was not longer visible in the granulation. Next, an initially homogeneous chromosphere with 130 grid layers and a vertical extent of $1.31\cdot 10^{9}$ m was added on top. The chromosphere evolves much faster than the layers below and is therefore only added in the last step in order to save computational time. The final model consists of $280\times 280\times 280$ grid cells and has a total spatial extent of $9.2\cdot 10^{9}\,\mathrm{m}\,\times\,9.2\cdot 10^{9}\,\mathrm{m}\,\times\,4.1\cdot 10^{9}\,\mathrm{m}$. The grid spacing is constant in horizontal direction ($\Delta x=\Delta y=3.29\cdot 10^{7}$ m) while the vertical grid cell size decreases from $\Delta z=5.73\cdot 10^{7}$ m) at the lower boundary in the convection zone to $\Delta z=9.66\cdot 10^{6}$ m in the atmosphere, i.e. for all heights above $z=0$ m (except for a slight increase in the topmost layers). The simulation with a chromosphere is then advanced for $5.7\cdot 10^{5}$ sec (157 hours) of stellar time. The model maintains an effective temperature of $T_{\mathrm{eff}}\approx(4010\pm 13)$ K. For comparison, we also use a model without a chromosphere, i.e., the relaxed model obtained before adding the extra 130 grid layers. Note that the entropy flux at the bottom boundary of the chromospheric model was adjusted to obtain the same total emergent radiative flux as that of the model without the chromosphere. In effect, this ensured that the effective temperatures of the two models are nearly identical.

2.3 Intensity synthesis

The radiative transfer code LINFOR3D111http://www.aip.de/Members/msteffen/linfor3d is used to calculate the emergent continuum intensity at different wavelengths ranging from the ultraviolet to the millimeter range (see Sect. 3.2 for the results). LINFOR3D is originally based on the Kiel code LINFOR/LINLTE and solves the radiative transfer equation in detail in 3D for an input atmosphere under the assumption of local thermodynamic equilibrium (LTE).

The same model was used as input for the MULTI_3D code, i.e. the column-by-column version of the original MULTI_3D code (Carlsson 1986), which provides the detailed solution of the radiative transfer equation in non-LTE (non-local thermodynamic equilibrium). The output contains intensity cubes ($I=f(x,y,\lambda)$) for several spectral lines of hydrogen (e.g., H$\alpha$ and Lyman $\alpha$) and singly ionized calcium (e.g., Ca II H, K, and the infrared triplet) for each selected snapshot.

2.4 Computation of the spectral energy distributions

The spectral energy distribution (SED) of the emergent radiation field is calculated using the NLTE3D code (Steffen et al. 2015). The code is designed to compute departure coefficients for the different energy levels of a given model atom222The departure coefficient, $b_{\rm i}$, of atomic level $i$ is defined as the ratio of population number densities, $b_{\rm i}=n_{\rm i}^{\rm NLTE}/n_{\rm i}^{\rm LTE}$, obtained under the assumptions of non-local thermodynamic equilibrium, NLTE, and local thermodynamic equilibrium, LTE.. For this purpose, NLTE3D calculates mean intensities at each grid point, which we then use to compute the SED for the given model atmosphere, i.e., the radiative flux at its top, at different wavelengths.

We compute SEDs for two model atmospheres, one with the chromosphere and one without it (see Sect. 2.2). Note that in both cases the two models shared identical atmospheric parameters, chemical composition, opacities, and equation of state. The SEDs were computed in the $150-10000$ nm wavelength range. To account for the line opacities, we use LITTLE opacity distribution functions (ODFs) from the ATLAS9 model atmosphere package (Castelli et al. 2003). Continuum opacities are taken into account by using the CO5BOLD IONOPA package (see Steffen et al. 2015, for details). The SEDs computed using the two model atmospheres are shown in Fig. 5.

3 Results

The initially homogeneous atmosphere evolves quickly towards a new dynamic equilibrium state. The first shock waves have propagated through the whole chromosphere and reach the upper boundary of the model after 210 000 s. From then on, the whole chromosphere is characterised by a shock-induced dynamic pattern, which is visible in many different quantities such as, e.g., gas temperature and velocity (Fig. 2). The resulting chromosphere exhibits many features known from earlier (non-magnetic) solar models (e.g., Wedemeyer et al. 2004; Skartlien et al. 2000) although the details are quite different. Propagating shock waves act as a structuring agent and produce a hot mesh of filaments, which can be seen in horizontal and vertical cross-sections through the selected model snapshot presented in Fig. 1 and Fig. 2, respectively. The resulting complex and dynamic topology exhibits variations on a large range of spatial scales, which are in general larger than in solar models but also involve scales as small as the extent of the intergranular lanes in the photosphere below. The thermal structure of the atmosphere and the emergent continuum intensity are addressed in Sect. 3.1 and Sect. 3.2, respectively.

3.1 Atmospheric structure and dynamics

The horizontal cross-sections for a selected model snapshot after 146.8 hours simulated time in Fig. 1 demonstrate how the atmospheric structure changes with height in terms of gas temperature, logarithmic gas pressure, and velocity. The convection zone (see the panels a-c) is characterised by mostly hot upwelling gas with a mesh of cool sinking gas. The cool gas is mostly concentrated in plumes that have their origin in intergranular vertices in the photosphere above. These plumes can also be seen as areas of lower gas pressure and in the velocity, which sometimes reaches downward speeds in excess of –15 km s*-1*. These flows can thus become supersonic and reach a Mach number of up to $\sim 2$. The average gas temperature at the bottom of the models is $\sim 17,500$ K and decreases monotonically to the effective temperature at the height where $\log\tau_{\mathrm{Ross}}=0$, i.e., at the transition to the photosphere.

The low photosphere (see the panels d-f in Fig. 1) is characterised by a granulation pattern with hot granules, where gas is rising to the surface, and narrow intergranular lanes, where the cooled gas sinks down again into the convection zone. The typical size of a granule is around $1.3\,-\,1.8\cdot 10^{9}$ m so that about 5 - 7 of such granules fit next to each other in the computational box. However, the distribution of granule sizes in this layer seems to be otherwise continuous with larger and smaller granules occurring, too. The intergranular lanes typically have widths of $0.2\,-\,0.3\cdot 10^{9}$ m with larger downdraft areas located at granule vertices. There, the conservation of momentum in the downflowing plasma results in vortex flows, which seem to be an integral part of stellar surface convection (see Wedemeyer-Böhm et al. 2012 and Wedemeyer et al. 2013, for vortex flows occurring in models of the Sun and M-dwarfs, respectively).

Shock formation.

The upwards propagating wave fronts steepen into shocks above $4\cdot 10^{8}$ m and sometimes as low as $3\cdot 10^{8}$ m and reach peak temperatures on the order of up to 5 000 K. The gas in the wake of shock waves is expanding adiabatically, which results in reduced gas temperatures as low as 2 000 K. The average gas temperature in the chromosphere is about 2 500 K. The velocities in shock fronts reach typical values on the order of 10 km s*-1*, with more extreme values of up to more than $20$ km s*-1*. As can be seen in Fig. 2f, the corresponding Mach numbers are typically on the order of 3 to 4 but can also reach values of 6 and - in very extreme, localised cases - up to 9. Upward propagating wave fronts can clearly be seen in vertical cross-section in Fig. 2e as regions of high upward directed velocity. It is also imminent that these wave fronts usually run into material that is falling down from above at high speeds, which often exceed $-20$ km s*-1*. The situation is also illustrated for a selected column (i.e., a single horizontal position in the computational model box) in Fig. 3. The column is taken from the vertical $x-z$ cut shown in Fig. 2 at a horizontal position of $x=3.3\cdot 10^{9}$ m. In that column, two major shock fronts can be seen at heights of approximately of $z=0.48\cdot 10^{9}$ m and $z=1.30\cdot 10^{9}$ m. The shock fronts can clearly be seen as temperature spikes and characteristic sawtooth profiles in the vertical velocity. The vertical velocity and the corresponding Mach number in Fig. 2c illustrate the existence of supersonic downflowing material. The fast downflows are caused by material that has been transported up by previous shock waves and then falls down again under the influence of gravity. These downflows themselves contribute to the continued formation of shocks throughout the whole model chromosphere. The same process has been found in earlier solar models (e.g., Carlsson & Stein 1995; Wedemeyer et al. 2004).

Gas temperature distribution.

An important consequence of the dynamic shock pattern in the chromosphere is that the thermal and kinetic state of the gas cannot be described well with average values only. In other words, the complicated and intermittent structure of the modelled atmosphere cannot be described correctly with a simple one-dimensional and static model atmosphere. The same is also true for the modelled convection zone. This result becomes obvious when looking at the histograms in Fig. 4. The displayed distributions for the same model snapshot as in Figs. 1 and 2 are very similar to the distributions for time steps spanning the second half of the chromospheric simulation. Hereafter, we refer to the latter. The gas temperature distribution for the chromosphere (Fig. 4g) has a peak at $T_{\mathrm{gas}}\approx 2\,450$ K whereas the average gas temperature is 2 670 K. The true peak of the distribution is still within one standard deviation (340 K) from the average value but there is a tail stretching to the maximum value of 4 600 K, which is not well presented by the average value at all. The high-temperature tail is due to the narrow chromospheric shock fronts, whereas the pronounced peak of the distribution is due to the cooler post-shock regions with temperatures down to $\sim 2\,000$ K. Due to the non-linear dependence of gas temperature and emergent intensity, the high-temperature tail significantly affects the average temperature derived from observed intensities as we will discuss in Sect. 4.

Velocity distribution.

The velocity distributions for selected heights in the model are shown in the rightmost column of Fig. 4. Again, the distributions for the selected model snapshot (red) are very close to the distributions for the whole second half of the simulation sequence, indicating that the selected snapshot is representative. For all selected heights, each distribution has a pronounced peak at downward directed velocity ($v_{z}<0$) and a peak for upwards directed velocity ($v_{z}>0$). In the convection zone (see Fig. 4c), roughly as many grid cells with downward velocities as grid cells with upward velocities are found, resulting in an average velocity close to zero. However, the peak for upward velocities is narrower with a maximum around $<|v|\times v_{z}/|v_{z}|>\approx 2$ km s*-1* whereas the downward velocities span a broader range of values with a maximum close to $<|v|>\approx-3$ km s*-1* and a tail reaching beyond $-10$ km s*-1*. The distribution is consistent with the rather smoothly upwelling hot gas and the cool gas shooting back into the convection zone in plumes. The same is in principle true for the “surface” layer at optical depth unity ($\log\tau_{\mathrm{Ross}}=0.0$, see Fig. 4f) although the distribution is more symmetric in the sense that the downward velocity peak is similar (but not identical) to the upward velocity part. The downward velocity distribution peaks around $<|v|>\approx-5...-6$ km s*-1* and has a tail extending well beyond $-10$ km s*-1*, even reaching $-15$ km s*-1* in extreme cases. The upward velocity distribution is slightly narrower, peaks around $<|v|>\approx+3...+4$ km s*-1* and has a tail extending well beyond $+10$ km s*-1*.

The situation is reversed in the chromosphere ($\log\tau_{\mathrm{Ross}}=-6.0$, see Fig. 4i) compared to the convection zone. In the chromosphere, slightly more grid cells exhibit downward velocities due to material falling down after having been lifted up by previous shock wave trains. The downward velocity part of the distribution peaks around $<|v|>\approx-10...-7$ km s*-1* with extreme values reaching $-20$ km s*-1*. Consequently, the arithmetic velocity average is negative, namely $<|v|>\approx-2$ km s*-1*. The upward velocity part of the distribution, which is due to upwards propagating shock waves, peaks around $<|v|>\approx+7...+10$ km s*-1* with extreme values reaching $+20$ km s*-1*. In summary, the model chromosphere is found to be highly dynamic with the matter often travelling at supersonic speeds and the whole atmospheric structure thus varying significantly on rather short time scales. At a typical speed of 10 km s*-1*, a shock wave (or equivalently downfalling material) propagates a distance of $10^{9}$ km, which is on the order of the thickness of the modelled chromospheric layer, in on the order of $10^{5}$ s.

3.2 Properties of the radiation field

Spectral energy distributions.

A comparison of SEDs computed in Sect. 2.4 and shown in Fig. 5 reveals that the red giant model with the chromosphere emits significantly more flux in the UV ($<270$ nm) than does the model without it. This is caused by the increasing average emissivity temperature towards the lower optical depths in the chromosphere. Obviously, the chromosphere makes an important contribution towards the emergent flux in the UV so that the model without a chromosphere naturally produces a significantly lower UV flux than the model with a chromosphere. This expected finding illustrates the need for detailed models that incorporate all atmospheric layers.

Continuum intensity maps.

In Figure 6 continuum intensity maps at different wavelengths for the same model snapshot as in Figs. 1 and 2 are presented. All maps represent a synthetic observation at disk center ($\mu=1.0$), i.e., as if seen spatially resolved from top. The maps for the wavelengths 300 nm, 500 nm, 800 nm, 1.6 $\mu$m all clearly exhibit the granulation pattern in the model although mapping slightly different height ranges around $\log\tau_{\mathrm{Ross}}=0$. Consequently, the contrasts of the maps decrease from 53.85 % at $\lambda=300$ nm to only 10.93 % at $\lambda=1.6\,\mu$m. The average intensity is shown as function of wavelengths in Fig. 6f. It has a peak at 500 nm.

The continuum intensity at millimeter wavelengths is another promising way to map the chromospheric layers. The large diagnostic potential of the Atacama Large Millimeter/submillimeter Array (ALMA) for this purpose has been demonstrated for the Sun (see Wedemeyer et al. 2016, and references therein) and other stars, e.g., $\alpha$ Cen (Liseau et al. 2015) and Mira (Vlemmings et al. 2015). In Fig. 6e a corresponding intensity map for a wavelength of 1.0 mm is shown. In contrast to the other maps at shorter wavelengths, the continuum intensity at 1.0 mm emerges from the chromosphere. The map therefore exhibits a filamentary pattern with apparent cells with diameters on the order of $1-2\cdot 10^{9}$ m and filament widths on the order of $0.1\cdot 10^{9}$ m. The pattern in intensity corresponds to the chromospheric shock fronts and cooler post-shock regions.

Spectral lines.

In Figure 7, the consequences of the intermittent nature of the model chromosphere are illustrated for a few commonly used spectral lines, namely Ca II K, the Ca II infrared triplet line at $\lambda=854$ nm, and H$\alpha$. The figure shows spatially averaged spectral line profiles and corresponding value ranges next to intensity maps for the line cores. Since the line cores are expected to form highest in the model chromosphere, the resulting line core maps should exhibit the intermittent structure already seen directly in the maps for gas temperature, pressure and velocity (see Fig. 1g-i). However, the different spectral lines are sensitive to different aspects of the model chromosphere due to non-linear dependencies in the line formation processes. This effect is particular obvious when comparing the Ca II line core maps in Fig. 7c and f with the continuum intensity map at $\lambda=1.0$ mm (Fig. 6). The shock-induced pattern of bright filamentary threads and dark regions is as clearly seen in the mm map as in the gas temperature map. The Ca II K map emphasizes more the hottest filamentary parts, whereas the picture is less clear in the Ca II 854 nm line core map. The contrast, i.e., the standard deviation of the quantity divided by its average, for the chromospheric gas temperature as shown in Fig. 1g is 14.9 %, whereas the contrasts for the Ca II K line and the Ca II 854 nm line core intensity maps are 274.6 % and 64.5 %, respectively. The mm radiation provides a more direct measure of the actual gas temperatures in the model chromosphere, which can be calculated under the assumption of local thermodynamic equilibrium (LTE), whereas the Ca II line cores are subject to non-LTE effects and a more complicated relationship between the state of the atmospheric gas and the line core intensity.

The spectral line profiles therefore vary significantly for the different locations. The intensity range covered by the individual Ca II K line profiles is illustrated in Fig. 7. For some horizontal positions, the line profiles exhibit Ca II K2 peaks in the core and equivalent for the Ca II 854 nm line core as can be seen from the upper 99 % percentile (upper dashed lines in panels b and e, respectively). However, such features are not as clearly visible for many other positions so that, despite the large spatial variations, the spatially averaged line profile does not exhibit a strong central emission reversal peak. While this makes it difficult to deduce the existence of shock waves from averaged spectra, as they would be observed, it obviously does not mean that shock waves do not exist. Fitting an observed spectrum with a one-dimensional model atmosphere would in this case lead to wrong conclusions by underestimating the extent to which shock waves occur.

The last row in Fig. 7 shows the synthetic intensities for the H$\alpha$ for the same model snapshot as discussed above. This spectral line is commonly used as chromospheric diagnostics and as activity indicator. Observations in H$\alpha$ on the Sun typically exhibit a complicated topology with fibrils that essentially outline the magnetic field (see, e.g., De Pontieu et al. 2007). The lack of such fibrils in the H$\alpha$ line core map for the red giant model shown in Fig. 7i is not unexpected because no magnetic fields are included in the model. Instead, the resulting intensity maps are dominated by photospheric contributions and thus mostly show the granulation pattern, which can be seen in the continuum intensity images in Fig. 6, too. The information that can be retrieved from analysing the H$\alpha$ line for this particular non-magnetic model is therefore quite different from the usual diagnostic use of the H$\alpha$ line.

4 Discussion

4.1 Preliminary comparison to observations

In Figure 7, the synthetic spectra for the presented red giant model are compared to observations of Aldebaran ($\alpha$ Tau), namely HARPS data333The HARPS data were obtained from http://www.blancocuaresma.com/s/. (Blanco-Cuaresma et al. 2014), UVES and VUES (Dobrovolskas et al. in prep.) data. A pipeline-reduced UVES spectrum (R$\sim$80 000) was taken from the UVES POP spectral library (Bagnulo et al. 2003). It should be noted that are remaining uncertainties regarding the scale as compared to the synthetic spectrum due to the difficulties in precisely determining the continuum in the observed spectrum. The VUES spectrum was obtained at Moletai Astronomical Observatory (Lithuania) using the 1.65 m telescope and Vilnius University Echelle Spectrograph (VUES Jurgenson et al. 2016) in high-resolution mode (R$\sim$67 000), taking $4\times 60$ sec exposures. The spectrum was reducing using the standard VUES pipeline procedure, and was further continuum-normalized using a synthetic 1D LTE spectrum of Aldebaran. Apart from omission of the line blends in the synthetic spectra, there are substantial differences in terms of spectral line depths and widths for all three spectral lines, which are not unexpected for several reasons. First of all, the numerical model has an effective temperature of $T_{\mathrm{eff}}\approx(4010\pm 13)$ K, which is thus 83 K higher than the 3927 K for Aldebaran as stated by Blanco-Cuaresma et al. (2014) (cf. Blackwell et al. 1991). More importantly, the exploratory model is still lacking physical ingredients that would influence the shapes of spectral lines formed in the chromosphere. In particular, the numerical simulation is purely hydrodynamic and thus represents an artificially quiet state while magnetic fields have recently been detected in $\alpha$ Tau and similar other stars (Aurière et al. 2015). An ultimate next, although computationally expensive step, will be the inclusion of a magnetic field, non-equilibrium ionization and excitation of hydrogen and other species, which should lead to more realistic models and thus may better fit the observational data (see also Sect. 4.3). Despite this, one should also keep in mind that modeling of hydrogen line formation in red giant atmospheres is a notoriously difficult task and that current models are consequently still not capable of reproducing observational data satisfactorily (see, e.g., recent study of Bergemann et al. 2016, and discussion therein). It may therefore still take some time until sufficiently realistic models will become available. Furthermore, the model presented here only includes a photosphere and a chromosphere whereas an observed spectrum naturally contains possible contributions from a corona, circumstellar gas, incl. stellar wind (Harper et al. 2011), which in principle could (partially) account for differences between the modelled and observed spectra.

It is nevertheless interesting to see that already the simplified but dynamic model chromosphere exhibits downward velocities that often exceed the upward velocities. Consequently, the line cores are accordingly shifted for lines of sight with a strong downward velocity component, which also leaves an imprint in the averaged spectra in Fig. 7. This effect is known from the solar spectrum (cf. Peter & Judge 1999, and references therein). A small net downward velocity on the order of $1-2$ km/s has also been detected in the low chromosphere of $\alpha$ Tau (Robinson et al. 1998), which may hint at the existence of shock waves.

More detailed, quantitative comparisons with observations are essential for the further improvement of numerical models and thus ultimately for a correct interpretation of the observations. For this purpose, complementary diagnostics that probe different aspects of stellar atmospheres should be combined. That includes different spectral lines and continua that are formed in different atmospheric layers and that are sensitive to different properties of the atmospheric gas (e.g., density, temperature, magnetic field). The centre-to-limb variation as far as it can be derived from stellar observations is another important test. For instance, Richichi et al. (2017) used lunar occultations of $\alpha$ Tau to investigate the limb-darkening characteristics and found substantial asymmetries in the photospheric brightness profile. Such measurements could be directly compared to corresponding synthetic observables (e.g., limb darkening profiles) as it will be presented in forthcoming publications.

4.2 Average temperature stratification.

As mentioned above, the shock waves are clearly visible as sawtooth-like profiles in the gas temperature and velocity (Figs. 1-3). The shock wave signature, however, is intermittent and thus not visible in the horizontal arithmetic average of a single snapshot over even after averaging over a few time steps. The average gas temperature on the geometric height scale ($z$), as shown in Fig. 8a-b, appears to remain roughly constant throughout the model chromosphere. Like for the solar models, there is no distinct temperature minimum visible in the averaged temperature stratification. Also, averaging on a column mass scale (Fig. 8c-d) or optical depth scale (Fig. 8e-f) produces no substantial temperature rise. On the other hand, the shaded areas in all panels illustrate the large range of temperature values, which is most extreme when looking at the temperature as function of optical depth in Fig. 8f. The extended value ranges, in line with the histograms in Fig. 4g, strongly suggest that the chromospheric gas temperature distribution cannot be described sufficiently by means of a simple arithmetic average alone. In the case of an atmosphere with significant spatial and/or temporal temperature variations, deriving an average temperature stratification by adjusting a model atmosphere to reproduce an intrinsically spatially average stellar spectrum only bears limited if not potentially misleading information. While this statement was proven to be true for the Sun (as described below), it should be noted that the model presented here is still of exploratory nature and more advanced simulations in the future are needed for a more realistic assessment of the thermal structure of red giant chromospheres.

The non-linear dependence of the source function and thus the emergent intensity on the gas temperature, as it is for instance the case in the ultraviolet, shifts the derived average gas temperature towards high temperatures in shock fronts, which enter the average with a higher weight than the cold temperatures in the cooler post-shock regions. The average emissivity temperature (cf. Wedemeyer et al. 2004)

[TABLE]

takes this effect into account, although rather crudely. The average emissivity temperature is shown as red dot-dashed lines in all panels of Fig. 8. In case of the Sun, this effect explains apparently contradicting temperature stratifications derived from UV continua and lines on the one hand and infrared carbon monoxide lines on the other hand (see, e.g., Ayres & Testerman 1981; Vernazza et al. 1981; Carlsson & Stein 1997, and references therein). The high-temperature component of the chromospheric temperature distribution is less pronounced in the red giant model presented here as compared to (non-magnetic) solar models (e.g., Wedemeyer et al. 2004) so that the influence of hot chromospheric gas on the average temperature is not as strong. Nevertheless, the emissivity temperature average produces a slightly higher average temperature in the chromosphere on the geometrical height scale (Fig. 8b), a notable chromospheric temperature rise on the column mass scale (Fig. 8d), and temperatures that are higher by up to 2 100 K on the optical depth scale (Fig. 8f). The resulting difference between the averages on the optical depth scale is particularly remarkable because the average emissivity temperature exceeds the 95 % percentile in the uppermost layers, indicating that a minority of extraordinarily hot locations in the model chromosphere strongly influence the resulting average due to the non-linear relation between the observable intensity and the actual gas temperature. It is this effect that has to be considered in detail when trying to derive meaningful atmospheric temperature stratification from intrinsically averaged stellar observations.

The model atmosphere derived by McMurry (1999) for $\alpha$ Tau is shown in Fig. 8 for comparison. The semi-empirical, plane-parallel, hydrostatic, one-component atmosphere model was constructed to reproduce observations of optical and ultraviolet photospheric and chromospheric spectral lines including, e.g., Ca II H and K, C I, C II, Si II, Mg II and C IV. There are small offsets of $\Delta z=7\,10^{7}$ m and $\Delta lgm_{c}=-0.25$, respectively, between the 3D numerical model and the McMurry model, which is not unexpected since the height of optical depth unity and the origin of the column mass scale vary for different horizontal positions in the 3D model so that the average scales depends on the exact way the average is calculated. On the column mass scale, the arithmetically averaged temperature from our simulation agrees well with the McMurry model for values of $\log m_{\mathrm{c}}>-1$. For lower values, i.e., above the low chromosphere the McMurry temperature stratification features a sharp chromospheric temperature rise and temperatures exceeding those in our model for $\log m_{\mathrm{c}}<-2$ and $z>1.4\,10^{9}$ m, respectively. The McMurry model also includes a transition region at $z=6\,10^{9}$ m with a jump to $10^{5}$ K which lies well above the upper boundary of the 3D simulation presented here. Although the 3D model does not exhibit a similarly strong temperature rise in the chromosphere, at least the average emissivity temperature begins to rise roughly in the same layer as the McMurry model although not as pronounced. As discussed above, plotting the temperatures on the optical depth scale444The optical depth scale for the McMurry model was calculated by converting the column mass scale with the help of the relation between optical depth and column mass in the 3D model. The resulting scale should be reasonably accurate for preliminary comparisons as presented here. results in a stronger apparent chromospheric increase of the average emissivity temperature. The corresponding 95 % percentile on the optical depth scale gives the impression of a sharp temperature rise similar to the McMurry model although this behaviour only applies the hottest chromospheric grid cells in the 3D model. Based on the exploratory 3D model, we conclude that (i) the chromospheric temperature rise has to be interpreted in view of potential averaging effects for a spatially intermittent atmosphere and that (ii) the discrepancy in temperature in the upper layers implies that the 3D simulations are still lacking physical mechanisms that are important for atmospheric heating. Overall the situation resembles the state of development of numerical model atmospheres of the Sun in the past for which the apparent controversy regarding the chromospheric temperature structure – in particular regarding its spatially and temporally intermittent nature – has been largely resolved by now (see, e.g., Avrett & Loeser 2008; Gudiksen et al. 2011; Rutten 2007; Vecchio et al. 2009; Wedemeyer-Böhm et al. 2009, and references therein).

4.3 Magnetic fields

Magnetic fields have been omitted for the first simulation presented here but will be considered for future models. Based on experience with CO5BOLD models with chromospheres for the Sun and other stellar types (Freytag et al. 2012; Wedemeyer et al. 2013), we expect that the inclusion of magnetic fields will change the chromospheric structure seen in the models significantly. While shock waves will also be excited in magnetohydrodynamic models, the details of their interaction with magnetic fields in the chromosphere will depend on the magnetic field strength and field topology in the model and change gradually from the non-magnetic shock-induced pattern like in the model presented here towards a magnetically dominated chromosphere with a mixture of open magnetic field lines and closed magnetic loops. Effects like, e.g, wave guiding, wave mode conversion, and initiation of torsional and rotational motions are expected and will not only affect the atmospheric structure but also the energy transport and thus the overall activity level of the model. In general, it is thought that stellar chromospheres are produced by a mixture of acoustic and magnetic heating mechanisms (e.g., Narain & Ulmschneider 1996). Although the results of the red giant simulation presented in this study reveal noticeable heating of its chromospheric layers due to shock waves, the importance of magnetic heating still needs to be assessed and will be considered in our future work. Interestingly, although we detect vigorous shock wave activity in the chromospheric layers of the red giant modeled here, red giants that are thought to be chromospherically inactive do not show spectroscopic signatures of shock waves (e.g., Judge & Carpenter 1998). This finding may suggest a magnetic origin of the heating of red giant chromospheres but, on the other hand, could simply mean that the shock signatures are by nature difficult to detect due to the intrinsic averaging over the whole stellar disk as illustrated for the Ca II lines in Fig. 7. On the other hand, there is observational evidence that suggests the presence of magnetic fields in the atmosphere of $\alpha$ Tau (e.g., Harper et al. 2011; Aurière et al. 2015, and references therein), which thus motivates the development of magnetohydrodynamic 3D simulations of this star.

5 Conclusions

The three-dimensional hydrodynamic model presented here is a first step towards more realistic model chromospheres for red giant stars that can support a meaningful in-depth analysis of stellar observations. Future models will include magnetic fields and eventually further physical ingredients like time-dependent hydrogen ionisation, which are important for the chromospheric gas properties and, e.g., the formation of a transition region (if existing for the modelled star) as natural upper boundary of a chromosphere. Despite the simplifications, the first model presented here already exhibits a very intermittent and dynamic chromosphere, which is in line with chromospheres modelled for other stellar types. The presented synthetic intensity maps for spectral lines and continua formed in the chromosphere clearly show that the exhibited spatial and temporal variations cannot be correctly described with a static one-dimensional model atmosphere. While 1D models have become very elaborate (see, e.g., Dupree et al. 2016), the continued development of adequate 3D (magneto-)hydrodynamical model atmospheres is ultimately needed.

Acknowledgements.

We thank V. Dobrovolskas for providing the VUES spectrum of Aldebaran prior to its publication. This work was supported by grant from the Research Council of Lithuania (MIP-089/2015). SW acknowledges support by grant from the Research Council of Lithuania (VIZ-TYR-158).