EDP Sciences
Highlight
Free Access
Issue
A&A
Volume 606, October 2017
Article Number A26
Number of page(s) 12
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/201730405
Published online 03 October 2017

© ESO, 2017

1. Introduction

The existence of chromospheres around red giant stars is suggested by a variety of observational evidence. Mostly, this evidence comes from indirect indicators related to stellar activity, 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α 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 to the presence of chromospheres (e.g. Dehaes et al. 2011; Harper et al. 2013). From a theoretical point of view, 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 and fluxes in the milimeter wavelength range (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 α Tau (McMurry et al. 1999). Detailed multi-component models with a realistic description of time-dependent shock waves require the use 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 (e.g. shock wave activity) that shape their structures.

In this work, we present the results of our first exploratory simulation of a 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, α 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 a rough approximation of the conditions expected for very quiet regions in the low chromosphere. While this is already of interest, α 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 the form of a look-up table, which is computed in advance under the assumption of thermodynamic equilibrium. For this table, partial ionisation of H and He, the formation and dissociation of H2, 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 ionisation degree in chromospheres have not been considered 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 factors are taken into account. A long characteristic 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 were 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 in several steps (see Sect. 2.1). The initial model atmosphere is based on a snapshot of a relaxed simulation that was computed using the following atmospheric parameters: Teff = 4000 K, log g = 1.5, and [ 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 × 109 m × 4.6 × 109 m × 2.8 × 109 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 × 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 quickly as possible. The sequence was run for one million seconds of stellar simulation time until the initial symmetry was no longer visible in the granulation. Next, an initially homogeneous chromosphere with 130 grid layers and a vertical extent of 1.31 × 109 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 × 280 × 280 grid cells and has a total spatial extent of 9.2 × 109 m × 9.2 × 109 m × 4.1 × 109 m. The grid spacing is constant in horizontal direction (Δx = Δy = 3.29 × 107 m), while the vertical grid cell size decreases from Δz = 5.73 × 107 m at the lower boundary in the convection zone to Δz = 9.66 × 106 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 × 105 s (157 h) of stellar time. The model maintains an effective temperature of Teff ≈ (4010 ± 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. We 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. This ensured that the effective temperatures of the two models are nearly identical.

2.3. Intensity synthesis

The radiative transfer code LINFOR3D1 was used to calculate the emergent continuum intensity at different wavelengths ranging from the ultraviolet to the millimetre range (see Sect. 3.2 for the results). LINFOR3D is 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 code (Carlsson 1986), which provides the detailed solution of the radiative transfer equation in non-LTE. The output contains intensity cubes (I = f(x,y,λ)) for several spectral lines of hydrogen (e.g. Hα and Lyman α) and singly ionised 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 atom2. 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). We 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 − 10 000 nm wavelength range. To account for the line opacities, we used the LITTLE opacity distribution functions (ODFs) from the ATLAS9 model atmosphere package (Castelli et al. 2003). Continuum opacities are taken into account by using the CO5BOLDIONOPA package (see Steffen et al. 2015, for details). The SEDs computed using the two model atmospheres are shown in Fig. 5.

thumbnail Fig. 1

Selected physical quantities in horizontal cross-sections through a single model snapshot (i.e. 3D model structure computed at a single instant in time) of the red giant model atmosphere after 146.8 h simulated time. The different columns show maps for the gas temperature (left), logarithm of gas pressure (middle), and the product of the absolute velocity and the sign of the vertical velocity component (right). Negative velocities are directed downwards/inwards (i.e. towards the stellar core) and positive velocities upwards/outwards. Each of these quantities is presented on planes of constant optical depth in the convection zone (top, log τRoss = 6.0), the photosphere (middle, log τRoss = 0.0), and the chromosphere (bottom, log τRoss = − 6.0). The streamlines in the right column follow the horizontal velocity field in the depicted plane.

Open with DEXTER

thumbnail Fig. 2

Vertical cross-sections along the x-axis at y = 0.577 × 109 m through the same model snapshot as in Fig. 1, displaying the following physical quantities: a) gas temperature, b) logarithmic gas pressure, c) logarithmic gas density, d) the difference between the local logarithmic optical depth and its horizontal average (log τRoss(x,z) − < log τRoss>x,y), e) absolute velocity multiplied by the sign of the vertical velocity component, and f) the Mach number. In panels e and d, streamlines for the flow field in the view plane and a dot-dashed contour for v = 0 km s-1, which thus divide upflows (v> 0) from downflows (v< 0), are plotted. The height of optical depth unity is represented by the black solid line around z = 0 m, which thus marks the boundary between photosphere and convection zone.

Open with DEXTER

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 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 the horizontal and vertical cross-sections through the selected model snapshot presented in Figs. 1 and 2, respectively. The resulting complex and dynamic topology exhibits variations on a wide 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 Sects. 3.1 and 3.2, respectively.

3.1. Atmospheric structure and dynamics

The horizontal cross-sections for a selected model snapshot after 146.8 hours simulated time (Fig. 1) demonstrate how the atmospheric structure changes with height in terms of gas temperature, logarithmic gas pressure, and velocity. The convection zone (see panels ac) 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 ~2. The average gas temperature at the bottom of the models is ~17 500 K and decreases monotonically to the effective temperature at the height where log τRoss = 0, i.e. at the transition to the photosphere.

The low photosphere (see panels df 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 × 109 m so that about 57 granules fit next to each other in the computational box. However, the distribution of granule sizes in this layer seems to be continuous with larger and smaller granules occurring as well. The intergranular lanes typically have widths of 0.2 − 0.3 × 109 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 upward propagating wave fronts steepen into shocks above 4 × 108 m and sometimes as low as 3 × 108 m and reach peak temperatures of up to 5000 K. The gas in the wake of shock waves expands adiabatically, which results in reduced gas temperatures as low as 2000 K. The average gas temperature in the chromosphere is about 2500 K. The velocities in shock fronts reach typical values of 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 of 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 be clearly seen in vertical cross-section in Fig. 2e as regions of high upward directed velocity. These wave fronts usually run into material that is falling down from above at high speeds, often exceeding − 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 xz cut shown in Fig. 2 at a horizontal position of x = 3.3 × 109 m. In that column, two major shock fronts can be seen at heights of approximately z = 0.48 × 109 m and z = 1.30 × 109 m. The shock fronts can be clearly 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 presence 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 1D 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 distributions shown 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 Tgas ≈ 2450 K, whereas the average gas temperature is 2670 K. The true peak of the distribution is still within one standard deviation (340 K) of the average value, but there is a tail stretching to the maximum value of 4600 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 ~2000 K. Owing to the non-linear dependence of gas temperature and emergent intensity, the high-temperature tail significantly affects the average temperature derived from observed intensities (see 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 (vz< 0) and a peak for upward directed velocity (vz> 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 | × vz/ | vz | ⟩ ≈ 2 km s-1, whereas the downward velocities span a broader range of values with a maximum close to ⟨ | v | ⟩ ≈ − 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. In principle, the same is true for the surface layer at optical depth unity (log τ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 | ⟩ ≈ − 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 | ⟩ ≈ + 3... + 4 km s-1, and has a tail extending well beyond + 10 km s-1.

The situation is reversed in the chromosphere (log τ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 | ⟩ ≈ − 10... − 7 km s-1 with extreme values reaching − 20 km s-1. Consequently, the arithmetic velocity average is negative, namely ⟨ | v | ⟩ ≈ − 2 km s-1. The upward velocity part of the distribution, which is due to upward propagating shock waves, peaks around ⟨ | v | ⟩ ≈ + 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 timescales. At a typical speed of 10 km s-1, a shock wave (or equivalently downfalling material) propagates a distance of 109 km, which is approximately the thickness of the modelled chromospheric layer, i.e. of the order of 105 s.

3.2. Properties of the radiation field

thumbnail Fig. 3

Vertical cross-section along the height axis z through the same model snapshot as in Figs. 1 and 2 at the horizontal position [ x,y ] = [ 3.31 × 109 m,0.58 × 109 m ]. The following physical quantities are displayed: a) gas temperature, b) vertical velocity, and c) Mach number. The dotted vertical lines mark the positions in two major shock fronts, where the shocks meet downfalling gas.

Open with DEXTER

thumbnail Fig. 4

Histograms for the horizontal cross-sections displayed in Fig. 1 for the same model snapshot (red) and for snapshots stretching over the second half of the simulation (grey/black). The columns are in the same order as in Fig. 1: gas temperature, logarithm gas pressure, and the product of the absolute velocity and the sign of the vertical velocity component, all on planes of equal optical depth. The rows are, from top to bottom: convection zone (log τRoss = 6.0), photosphere (log τRoss = 0.0), and chromosphere (log τRoss = − 6.0). The solid vertical lines mark the average values, whereas the dashed lines mark the average plus/minus one standard deviation.

Open with DEXTER

thumbnail Fig. 5

Top: spectral energy distributions (SEDs) of the red giant models with and without chromosphere. The blue line shows the average SED of the model without chromosphere computed using an ensemble of 20 model snapshots obtained at different instants in time. The red line shows the SED of the model with chromosphere computed using the single snapshot shown in Fig. 1. Bottom: relative difference between the flux of the models with and without chromosphere. The vertical scale is different in the wavelength region <270 nm.

Open with DEXTER

thumbnail Fig. 6

Emergent continuum intensity for a selected simulation snapshot at different wavelengths: a) 300 nm, b) 500 nm, c) 800 nm, d) 1.6 μm, and e) 1.0 mm. Panel f shows the average intensity (solid line with circles) and the corresponding intensity contrast (dot-dashed line with squares) of the maps displayed in panels a–e.

Open with DEXTER

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.

thumbnail Fig. 7

Results of radiative transfer calculations with the MULTI non-LTE code for the selected model snapshot for three different spectral lines: Ca II K (top row), Ca II 854.2 nm (middle row), and Hα (bottom row). For each row, the leftmost panel shows the horizontally averaged spectral line profile (thick solid line), the 10th and 90th percentiles (white dot-dashed with enclosed dark grey shaded area), the 1st and 99th percentiles (dashed with enclosed light grey shaded area), and the extreme values (dotted) found for that model snapshot. The same plots are repeated in the middle column, but for a narrower wavelength range around the line core. For comparison, normalised observed spectra for α Tau are shown as red lines in the left and middle column: UVES POP data (Bagnulo et al. 2003) in panels ab, VUES data in panels de (Dobrovolskas et al., in prep.), and HARPS data in panels gh (Blanco-Cuaresma et al. 2014). The red axis is to the right of the panels. The horizontally resolved intensity maps for the (nominal) line core wavelength are shown in the rightmost column for a view from the top of the model, thus corresponding to stellar disk centre.

Open with DEXTER

Continuum intensity maps.

In Fig. 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 centre (μ = 1.0), i.e. as if seen spatially resolved from top. The maps for the wavelengths 300 nm, 500 nm, 800 nm, and 1.6 μm all clearly exhibit the granulation pattern in the model although mapping slightly different height ranges around log τRoss = 0. Consequently, the contrasts of the maps decrease from 53.85% at λ = 300 nm to only 10.93% at λ = 1.6 μm. The average intensity is shown as function of wavelength in Fig. 6f. It has a peak at 500 nm.

The continuum intensity at millimetre 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. α 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 of the order of 1 − 2 × 109 m and filament widths of the order of 0.1 × 109 m. The pattern in intensity corresponds to the chromospheric shock fronts and cooler post-shock regions.

Spectral lines.

In Fig. 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 λ = 854 nm, and Hα. 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 Figs. 1gi). 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 particularly obvious when comparing the Ca II line core maps in Figs. 7c and f with the continuum intensity map at λ = 1.0 mm (Fig. 6). The shock-induced pattern of bright filamentary threads and dark regions is clearly seen in the mm map and equally well in the gas temperature map. The Ca II K map emphasises 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 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 99th percentile. However, these features are not as clearly visible for many other positions; 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 1D model atmosphere would in this case lead to the wrong conclusions by underestimating the extent to which shock waves occur.

The last row in Fig. 7 shows the synthetic intensities for Hα for the same model snapshot as discussed above. This spectral line is commonly used as a chromospheric diagnostic and as an activity indicator. Observations in Hα 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α 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 also be seen in the continuum intensity images in Fig. 6. The information that can be retrieved from an analysis of the Hα line for this particular non-magnetic model is therefore quite different from the usual diagnostic use of the Hα line.

4. Discussion

thumbnail Fig. 8

Average stratification of the gas temperature on the a), b) geometric height scale, c), d) column mass scale, and e), f) optical depth scale. The right column repeats the data in the left column, but only for the atmospheric layers (i.e. those located to the right of the vertical dashed line in the panels on the left). The arithmetic temperature averages are shown as solid black lines, whereas the corresponding average emissivity temperatures as defined in Eq. (1) are represented by the red dot-dashed lines. The light grey shaded areas mark the whole range of values between the minimum and maximum at each height or column mass or optical depth, whereas the dark grey areas represent the corresponding values between the 5th and 95th  percentiles, i.e. the majority of all values except for the most extreme ones. We note that these are distributions for the selected model snapshot only. For comparison, the α Tau model atmosphere by McMurry (1999) is plotted in all panels as blue dashed lines. Applying small offsets of Δz = 7 × 107 m and Δlgmc = − 0.25, respectively, produces a better match with our simulation and the shifted McMurry model (thick blue line). The optical depth scale for the McMurry is derived from matching the column mass scale to the optical depth scale in the 3D model.

Open with DEXTER

4.1. Preliminary comparison to observations

In Fig. 7, the synthetic spectra for the presented red giant model are compared to observations of Aldebaran (α Tau), namely HARPS data3 (Blanco-Cuaresma et al. 2014), and UVES and VUES data (Dobrovolskas et al., in prep.). A pipeline-reduced UVES spectrum (R ~ 80 000) was taken from the UVES POP spectral library (Bagnulo et al. 2003). It should be noted that there are 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 ~ 67 000), taking 4 × 60 s exposures. The spectrum was reduced using the standard VUES pipeline procedure, and was further continuum-normalised 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 Teff ≈ (4010 ± 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 α Tau and similar other stars (Aurière et al. 2015). The ultimate next, although computationally expensive step, will be the inclusion of a magnetic field, non-equilibrium ionisation, 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 modelling 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 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, and 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 the simplified but dynamic model chromosphere already 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 (Peter & Judge 1999, and references therein). A small net downward velocity of the order of 1 − 2 km s-1 has also been detected in the low chromosphere of α 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 α Tau to investigate the limb-darkening characteristics and found substantial asymmetries in the photospheric brightness profile. These measurements were able to be directly compared to corresponding synthetic observables (e.g. limb darkening profiles), as 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. 13). The shock wave signature, however, is intermittent and thus not visible in the horizontal arithmetic average of a single snapshot, even after averaging over a few time steps. The average gas temperature on the geometric height scale (z), as shown in Figs. 8a and b, appears to remain roughly constant throughout the model chromosphere. As is true for the solar models, there is no distinct temperature minimum visible in the averaged temperature stratification. Also, averaging on a column mass scale (Figs. 8c and d) or optical depth scale (Figs. 8e and f) produces no substantial temperature rise. On the other hand, the shaded areas in all panels illustrate the wide range of temperature values, which is most extreme when looking at the temperature as a 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 an 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, for instance as is 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 (Wedemeyer et al. 2004) (1)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 than 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 2100 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 95th 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 stratifications from intrinsically averaged stellar observations.

The model atmosphere derived by McMurry (1999) for α 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 Ca ii H and K, C i, C ii, Si ii, and Mg ii and C iv. There are small offsets of Δz = 7 × 107 m and Δlgmc = − 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 mc> − 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 mc< − 2 and z> 1.4 × 109 m, respectively. The McMurry model also includes a transition region at z = 6 × 109 m with a jump to 105 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, the average emissivity temperature begins to rise roughly in the same layer as the McMurry model although not as sharply. As discussed above, plotting the temperatures on the optical depth scale4 results in a stronger apparent chromospheric increase of the average emissivity temperature. The corresponding 95th 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 the chromospheric temperature rise has to be interpreted in view of potential averaging effects for a spatially intermittent atmosphere and that 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 significantly change the chromospheric structure seen in the models. 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 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 modelled 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 or, 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. However, there is observational evidence that suggests the presence of magnetic fields in the atmosphere of α 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 3D 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, for example, the formation of a transition region (if existing for the modelled star) as a 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 1D 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.


2

The departure coefficient, bi, of atomic level i is defined as the ratio of population number densities, , obtained under the assumptions of non-local thermodynamic equilibrium (NLTE) and local thermodynamic equilibrium (LTE).

3

The HARPS data were obtained from http://www.blancocuaresma.com/s/

4

The 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.

Acknowledgments

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

References

All Figures

thumbnail Fig. 1

Selected physical quantities in horizontal cross-sections through a single model snapshot (i.e. 3D model structure computed at a single instant in time) of the red giant model atmosphere after 146.8 h simulated time. The different columns show maps for the gas temperature (left), logarithm of gas pressure (middle), and the product of the absolute velocity and the sign of the vertical velocity component (right). Negative velocities are directed downwards/inwards (i.e. towards the stellar core) and positive velocities upwards/outwards. Each of these quantities is presented on planes of constant optical depth in the convection zone (top, log τRoss = 6.0), the photosphere (middle, log τRoss = 0.0), and the chromosphere (bottom, log τRoss = − 6.0). The streamlines in the right column follow the horizontal velocity field in the depicted plane.

Open with DEXTER
In the text
thumbnail Fig. 2

Vertical cross-sections along the x-axis at y = 0.577 × 109 m through the same model snapshot as in Fig. 1, displaying the following physical quantities: a) gas temperature, b) logarithmic gas pressure, c) logarithmic gas density, d) the difference between the local logarithmic optical depth and its horizontal average (log τRoss(x,z) − < log τRoss>x,y), e) absolute velocity multiplied by the sign of the vertical velocity component, and f) the Mach number. In panels e and d, streamlines for the flow field in the view plane and a dot-dashed contour for v = 0 km s-1, which thus divide upflows (v> 0) from downflows (v< 0), are plotted. The height of optical depth unity is represented by the black solid line around z = 0 m, which thus marks the boundary between photosphere and convection zone.

Open with DEXTER
In the text
thumbnail Fig. 3

Vertical cross-section along the height axis z through the same model snapshot as in Figs. 1 and 2 at the horizontal position [ x,y ] = [ 3.31 × 109 m,0.58 × 109 m ]. The following physical quantities are displayed: a) gas temperature, b) vertical velocity, and c) Mach number. The dotted vertical lines mark the positions in two major shock fronts, where the shocks meet downfalling gas.

Open with DEXTER
In the text
thumbnail Fig. 4

Histograms for the horizontal cross-sections displayed in Fig. 1 for the same model snapshot (red) and for snapshots stretching over the second half of the simulation (grey/black). The columns are in the same order as in Fig. 1: gas temperature, logarithm gas pressure, and the product of the absolute velocity and the sign of the vertical velocity component, all on planes of equal optical depth. The rows are, from top to bottom: convection zone (log τRoss = 6.0), photosphere (log τRoss = 0.0), and chromosphere (log τRoss = − 6.0). The solid vertical lines mark the average values, whereas the dashed lines mark the average plus/minus one standard deviation.

Open with DEXTER
In the text
thumbnail Fig. 5

Top: spectral energy distributions (SEDs) of the red giant models with and without chromosphere. The blue line shows the average SED of the model without chromosphere computed using an ensemble of 20 model snapshots obtained at different instants in time. The red line shows the SED of the model with chromosphere computed using the single snapshot shown in Fig. 1. Bottom: relative difference between the flux of the models with and without chromosphere. The vertical scale is different in the wavelength region <270 nm.

Open with DEXTER
In the text
thumbnail Fig. 6

Emergent continuum intensity for a selected simulation snapshot at different wavelengths: a) 300 nm, b) 500 nm, c) 800 nm, d) 1.6 μm, and e) 1.0 mm. Panel f shows the average intensity (solid line with circles) and the corresponding intensity contrast (dot-dashed line with squares) of the maps displayed in panels a–e.

Open with DEXTER
In the text
thumbnail Fig. 7

Results of radiative transfer calculations with the MULTI non-LTE code for the selected model snapshot for three different spectral lines: Ca II K (top row), Ca II 854.2 nm (middle row), and Hα (bottom row). For each row, the leftmost panel shows the horizontally averaged spectral line profile (thick solid line), the 10th and 90th percentiles (white dot-dashed with enclosed dark grey shaded area), the 1st and 99th percentiles (dashed with enclosed light grey shaded area), and the extreme values (dotted) found for that model snapshot. The same plots are repeated in the middle column, but for a narrower wavelength range around the line core. For comparison, normalised observed spectra for α Tau are shown as red lines in the left and middle column: UVES POP data (Bagnulo et al. 2003) in panels ab, VUES data in panels de (Dobrovolskas et al., in prep.), and HARPS data in panels gh (Blanco-Cuaresma et al. 2014). The red axis is to the right of the panels. The horizontally resolved intensity maps for the (nominal) line core wavelength are shown in the rightmost column for a view from the top of the model, thus corresponding to stellar disk centre.

Open with DEXTER
In the text
thumbnail Fig. 8

Average stratification of the gas temperature on the a), b) geometric height scale, c), d) column mass scale, and e), f) optical depth scale. The right column repeats the data in the left column, but only for the atmospheric layers (i.e. those located to the right of the vertical dashed line in the panels on the left). The arithmetic temperature averages are shown as solid black lines, whereas the corresponding average emissivity temperatures as defined in Eq. (1) are represented by the red dot-dashed lines. The light grey shaded areas mark the whole range of values between the minimum and maximum at each height or column mass or optical depth, whereas the dark grey areas represent the corresponding values between the 5th and 95th  percentiles, i.e. the majority of all values except for the most extreme ones. We note that these are distributions for the selected model snapshot only. For comparison, the α Tau model atmosphere by McMurry (1999) is plotted in all panels as blue dashed lines. Applying small offsets of Δz = 7 × 107 m and Δlgmc = − 0.25, respectively, produces a better match with our simulation and the shifted McMurry model (thick blue line). The optical depth scale for the McMurry is derived from matching the column mass scale to the optical depth scale in the 3D model.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.