© ESO, 2011

1. Introduction

The dynamical nature of the solar-surface layers, manifested for instance in granules and sunspots, has been known for a long time. With every improvement of ground-based or spaceborne instruments the complexity of the observed processes increased. Red supergiant (RSG) stars are among the largest stars in the universe and the brightest in the optical and near-infrared. They are massive stars with masses between roughly 10 and 25 M_⊙ with effective temperatures ranging from 3450 to 4100 K, luminosities of 20 000 to 300 000 L_⊙, and radii up to 1500 R_⊙ (Levesque et al. 2005). These stars exhibit variations in integrated brightness, surface features, and the depths, shapes, and Doppler shifts of spectral lines; as a consequence, stellar parameters and abundances are difficult to determine. Progress has been made using 1D hydrostatic models revising the T_eff-scale and reddening both at solar and Magellanic Clouds metallicities (Levesque et al. 2005, 2006; Massey et al. 2007; Levesque et al. 2007; Levesque 2010) but problems still remain, e.g. the blue-UV excess that may be caused by scattering by circumstellar dust or an insufficiency in the models, and the visual-infrared effective temperature mismatch (Levesque et al. 2006). Finally, RSGs eject massive amounts of mass back to the interstellar medium with an unidentified process that may be related to acoustic waves and radiation pressure on molecules (Josselin & Plez 2007), or to the dissipation of Alfvén waves from magnetic field, recently discovered on RSGs (Aurière et al. 2010; Grunhut et al. 2010), as early suggested by Hartmann & Avrett (1984); Pijpers & Hearn (1989); Cuntz (1997).

The dynamical convective pattern of RSGs is then crucial for the understanding of the physics of these stars that contribute extensively to the chemical enrichment of the Galaxy. There is a number of multiwavelength imaging examples of RSGs (e.g. α Ori) because of their high luminosity and large angular diameter. Several research teams (Buscher et al. 1990; Wilson et al. 1992; Tuthill et al. 1997; Wilson et al. 1997; Young et al. 2000) detected time-variable inhomogeneities on the surface of α Ori with WHT and COAST; Haubois et al. (2009) reported a reconstructed image in the H-band with two large spots; Ohnaka et al. (2009, 2011) detected possible convective motions in CO line formation layers, but it is not obvious whether this strong convective motion can be present in these layers (detected up to 1.3–1.4 stellar radii), where the convective energy flux is low; Kervella et al. (2009) resolved α Ori using diffraction-limited adaptive optics in the near-infrared and found an asymmetric envelope around the star with a bright plume extending in the southwestern region. Tuthill et al. (1997) reported bright spots on the surface of the supergiants α Her and α SCO using WHT, and Chiavassa et al. (2010c) on VX Sgr using VLTI/AMBER.

The effects of convection and non-radial waves can be represented by numerical multi-dimensional time-dependent radiation hydrodynamics (RHD) simulations with realistic input physics. Three-dimensional radiative hydrodynamics simulations are no longer restricted to the Sun (for a review on the Sun models see Nordlund et al. 2009) but cover a substantial portion of the Hertzsprung-Russell diagram (Ludwig et al. 2009). Moreover, they have been already extensively employed to study the effects of photospheric inhomogeneities and velocity fields on the formation of spectral lines and on interferometric observables in a number of cases, including the Sun, dwarfs and subgiants (e.g. Asplund et al. 1999; Asplund & García Pérez 2001; Asplund et al. 2009; Caffau et al. 2010; Behara et al. 2010; Sbordone et al. 2010), red giants (e.g. Chiavassa et al. 2010a; Collet et al. 2007, 2009; Wende et al. 2009), asymptotic giant branch stars (Chiavassa et al. 2010c), and RSG stars (Chiavassa et al. 2009, 2010b, 2011). In particular, the presence and the characterization of the size of convective cells on α Ori has been shown by Chiavassa et al. (2010b) by comparing an extensive set of interferometric observations ranging from the optical to the infrared wavelengths.

This paper is the fourth in a series aimed to present the numerical simulations of RSG stars with CO⁵BOLD and to introduce the new generation of RSG simulations with a more sophisticated opacity treatment. These simulations are very important for an accurate quantitative analysis of observed data.

2. 3D numerical simulations with CO⁵BOLD

2.1. Basic equations

The code solves the coupled equations of compressible hydrodynamics and non-local radiation transport\arraycolsep1.75ptThey describe the inviscid flow of density ρ, momentum ρv, and total energy e_ikg (including internal, kinetic, and gravitational potential energy). Further quantities are the velocity vector v, pressure P, and radiative energy flux . The latter is computed from the frequency-integrated intensity I for the gray treatment of opacity and for each wavelength group in the non-gray approach (Sect. 2.3).

The gravitational potential is spherical, (4)where M_pot is the mass of the star to be modeled (see Fig. 3, bottom right panel); r₀ and r₁ are free smoothing parameters: when , ; while when r → 0, ; and for r → ∞, . Typically, r₀ ≈ 0.2 R _⋆  and r₁ ≈ 1.2 R _⋆  for RSG simulations (Freytag et al. 2002).

2.2. The code

The numerical simulations described here are performed with CO⁵BOLD (COnservative COde for the COmputation of COmpressible COnvection in a BOx of L Dimensions, L = 2,3). It uses operator splitting (Strang 1968) to separate the various (explicit) operators: the hydrodynamics, the optional tensor viscosity, and the radiation transport.

The hydrodynamics module is based on a finite-volume approach and relies on directional splitting to reduce the 2D or 3D problem to one dimension. In the 1D steps an approximate Riemann solver of Roe-type (Roe 1986) is applied, modified to account for a realistic equation of state, a non-equidistant Cartesian grid, and the presence of source terms due to an external gravity field. In addition to the stabilizing mechanism inherent in an upwind-scheme with a monotonic reconstruction method (typically a piecewise-linear van Leer interpolation), a 2D or 3D tensor viscosity can be activated. This step eliminates certain errors of Godunov-type methods dealing with strong velocity fields aligned with the grid (Quirk 1994).

The equation of state uses pre-tabulated values as functions of density and internal energy (ρ,e_i → P,Γ₁,T,s). It accounts for Hi, Hii, H₂, Hei, Heii, Heiii and a representative metal for any prescribed chemical composition. The equation of state does not account for the ionization states of metals, but it uses only one neutral element to achieve the appropriate atomic weight (in the neutral case) for a given composition. Two different geometries can be used with CO⁵BOLD, which are characterized by different gravitational potentials, boundary conditions, and modules for the radiation transport:

• The box-in-a-star setup is used to model a statisticallyrepresentative volume of the stellar atmosphere with aconstant gravitation, where the lateral boundaries areperiodic, and the radiation transport module relies on aFeautrier scheme applied to a system of long rays (Freytaget al. 2002; Wedemeyeret al. 2004, and for analysis, e.g.,Caffau et al.2010).

• The star-in-a-box setup is used to model RSG stars of this work. The computational domain is a cubic grid equidistant in all directions, and the same open boundary condition is employed for all sides of the computational box.

Because the outer boundaries are usually either hit at some angle by an outgoing shockwave, or let material fall back (mostly with supersonic velocities), there is not much point in tuning the formulation for an optimum transmission of small-amplitude waves. Instead, a simple and stable prescription, that lets the shocks pass is sufficient. It is implemented by filling typically two layers of ghost cells where the velocity components and the internal energy are kept constant. The density is assumed to decrease exponentially in the ghost layers, with a scale height set to a controllable fraction of the local hydrostatic pressure scale height. The control parameter allows us to account for the fact that the turbulent pressure plays a significant role for the average pressure stratification. The acceleration through gravity is derived from Eq. (4). Within a radius r₀ the potential is smoothed (Eq. (4)). In this sphere a source term to the internal energy provides the stellar luminosity. Motions in the core are damped by a drag force to suppress dipolar oscillations. The hydrodynamics and the radiation transport scheme ignore the core completely and integrate right through it.

The code is parallelized with Open Multi-Processing (OpenMP) directives.

Fig. 1 Gray intensity on one side of the computational cube from the initial sequence of the model st35gm00n05 in Table 1. The axes are in solar radii. The artifacts caused by the mismatch between the spherical object and the Cartesian grid become less evident with time passing.

2.3. Radiation transport and opacities tables

The radiation transport step for RSGs’ simulations uses a short-characteristics method. To account for the short radiative time scale several (typically 6 to 9) radiative sub − steps are performed per global step. Each sub-step uses only three rays (e.g. in the directions (1,1,0), (1, −1,0), (0,0,1) or (1,0,1), (1,0, −1), (0,1,0) or (0,1,1), (0,1, −1), (1,0,0) or four rays (along the space diagonals (±1, ±1, ±1)). The different orientation sets are used cyclically. The radiation transport is solved for a given direction for both orientations. More irregular directions (and more rays) are possible but are avoided in the models presented here to save computational time. The radiation transport operator is constructed to be stable in the presence of strong jumps in opacity and temperature (even on a coarse grid).

The frequency dependance of the radiation field in the CO⁵BOLD models can be calculated with two approaches:

• the gray approximation, which completely ignores the fre-quency dependence, is justified only in the stellar interior andit is inaccurate in the optically thin layers. The Rosseland meanopacities are calculated as a function of pressure and temperature(T,P → κ_Ross) and are available in a 2D table. The necessary values are found by interpolation in a 2D table. It has been merged at around 12 000 K from high-temperature OPAL (Iglesias et al. 1992) data and low-temperature PHOENIX (Hauschildt et al. 1997) data by Hans-Günter Ludwig. The models that use these opacities are reported in Table 1.

• the more elaborate scheme accounting for non-gray effects, which is based on the idea of opacity binning (Nordlund 1982; Nordlund & Dravins 1990). The basic approximation is the so called multi-group scheme (Ludwig et al. 1994; Vögler et al. 2004). In this scheme, the frequencies that reach monochromatic optical depth unity within a certain depth range of the model atmosphere will be put into one frequency group. The opacity table used in this work is sorted to five wavelengths groups according to the run of the monochromatic optical depth in a corresponding MARCS (Gustafsson et al. 2008) 1D model. The corresponding logarithmic Rosseland optical depths are  +∞,0.0, −1.5, −3.0,  −4.5, −∞. In each group there is a smooth transition from a Rosseland average in the optically thick regime to a Planck average in the optically thin regime, except for the group representing the highest opacities, where the Rosseland average is used throughout. The implementation of non-gray opacities tables has been carried out by Hans-Günter Ludwig and we computed one non-gray model (Table 1).

3. Modeling RSG stars

The most important parameters that determine the type of the modeled star are

• the input luminosity into the core;

• the stellar mass that enters in Eq. (4) for the gravitational potential;

• the abundances that are used to create the tables for the equation-of-state and the opacities.

In addition, there are a number of parameters that influence the outcome of a simulation to some extent, i.e. the detailed formulation of the boundary conditions, the smoothing parameters of the potential, the numerical resolution, detailed settings of the hydrodynamics scheme and the additional tensor viscosity, choice of ray directions for the radiation transport, and the time-step control. The model presented in Freytag et al. (2002) relies on the same assumptions as the current ones. In the meantime, modifications of the code (e.g. restrictions to rays only along axes and diagonals avoiding some interpolations in the radiation transport module) and faster computers allow models with higher resolution, frequency-dependent opacity tables, and various stellar parameters. These new models show a significant increase in the number of small convective cells on the stellar surface.

The initial model is produced starting from a sphere in hydrostatic equilibrium with a weak velocity field inherited from a previous model with different stellar parameters (top left panel of Fig. 1). The temperature in the photosphere follows a gray relation. In the optically thick layers it is chosen to be adiabatic. The first frame of Fig. 1, taken just after one snapshot, displays the limb-darkened surface without any convective signature but with some regular patterns due to the coarse numerical grid and the correspondingly poor sampling of the sharp temperature jump at the bottom of the photosphere. The central spot, quite evident at the beginning of the simulation, vanishes completely when convection becomes strong. A regular pattern of small-scale convection cells develops initially and then, as the cells merge, the average structure becomes big and the regularity is lost. The intensity contrast grows. After several years, the state is relaxed and the pattern becomes completely irregular. All memory from the initial symmetry is lost. The influence of the cubical box and the Cartesian grid is strongest in the initial phase onset of convection. At later stages, there is no alignment of structures with the grid nor a tendency of the shape of the star to become cubical.

Fig. 2 Luminosity, temperature, and radius as a function of time for the simulations of Table 1: from left to right columns st35gm03n07, st35gm03n13, st36gm00n04, and st36gm00n05. The bottom panels are the ratio between turbulent pressure and gas pressure for different snapshots. The red vertical lines in all the panels is the approximative position of the radius from Table 1.

Fig. 3 Some spatially average quantities (i.e., spherical shells) for a snapshot of st35gm03n07 (Table 1). The red vertical line is the location of the radius. Top row, left panel: luminosities as function of radius in solar radii: continuous curve is the total luminosity; the dashed curve is the luminosity owing to pressure work; the dotted curve is the luminosity of kinetic energy; and the dash-dotted curve is the radiative luminosity. Top row, center: opacity. Top row, right: specific entropy. Central row, left panel: pressure scale height. Central row, center: characteristic radiative (continuous line) and hydrodynamical (dotted line) time scales. Central row, right: gravitational potential Φpot. Bottom row, left panel: random column opacity and radiative time scale for a given snapshot of the gray model st35gm03n07 (black) and non-gray model st35gm03n13 (red). The five colored curves correspond to the opacity groups. The radiative time scale is estimated for temperature fluctuations with spatial scales on the order of the grid size. Bottom row, center: characteristic radiative time scale for the gray (black) and non-gray (red) models. Bottom row, right: temperature profiles for the st36g00n04 (2553 grid points, black) and st36g00n05 (4013 grid points, red).

4. Model structure

Table 1 reports the simulations analyzed in this work. These models are the result of intensive calculations on parallel computers: the oldest simulation, st35gm03n07, used so far in the other papers of this series, took at least a year (this time may be difficult to estimate because of the queuing system and different machines used) on machines with 8–16 CPUs to provide nine years of simulated stellar time; more recently, we managed to compute the other simulations of Table 1 on a shorter time scale on a machine with 32 CPUs. The model st36gm00n04, which is slightly higher in numerical resolution, needed about two months to compute more than 20 stellar years, while st36gm00n05 (the highest numerical resolution so far) about six months for about three stellar years. Eventually, st35gm03n13 needed about eight months to compute ten years of stellar time despite its lower numerical resolution. The CPU time needed for non-gray runs does scale almost linearly with the number of wavelength groups.

4.1. Global quantities

Three-dimensional simulations start with an initial model that has an estimated radius, a certain envelope mass, a certain potential profile, and a prescribed luminosity as described above. During the run, the internal structure relaxes to something not to far away from the initial guess. We needed several trials to get the correct initial model for the latest set. The average final stellar parameters are determined once the simulation has ended (Fig. 2). For this purpose, we adopted the method reported in (Chiavassa et al. 2009, hereafter Paper I). The quantities are averages over spherical shells and over time, and the errors are one-sigma fluctuations with respect to the average over time. We computed the average temperature and luminosity over spherical shells, T(r), and L(r). Then, we searched for the radius R for which (5)where σ is the Stefan-Boltzmann constant. Eventually, the effective temperature is T_eff = T(R). As already pointed out in Paper I, the gray model st35gm03n07 (top row of Fig. 2) shows a drift in the first two years (−0.5% per year) before stabilisation in the last 2.5 years. This drift is also visible in the radius of the non-gray model st35gm03n13 (second row) even though it is weaker (less than 0.3%). The luminosity fluctuations are on the order of 1% and the temperature variations  <0.3% for both simulations. It must be noted that the last snapshot of the gray simulation was used as the initial snapshot for the non-gray model. The two hotter models, st36gm00n04 and st36gm00n05, are also somehow related because they both use a gray treatment of opacities and st36gm00n05 has been started using the last snapshot of st36gm00n04. The temperature (0.5%) and luminosity (2.2%) fluctuations of the higher resolution model are slightly higher than the temperature (0.4%) and luminosity (1.5%) fluctuations of st36gm00n04. None of the simulations shows a drift of the radius in the last years of evolution.

The bottom row in Fig. 2 shows the ratio between the turbulent pressure and the gas pressure (defined as in Sect. 4.2). This quantity shows that in the outer layers, just above the stellar radius, the turbulent pressure plays a significant role for the average pressure stratification and the radial velocities resulting from the vigorous convection are supersonic.

Figure 3 displays some quantities spatially averaged over spherical shells for a snapshot of the gray simulation st35gm03n07 in Table 1. Figure 4 shows three-dimensional views of some quantities for the same snapshot and simulation. Radiation is of primary importance for many aspects of convection and the envelope structure in a RSG. It does not only cool the surface to provide a somewhat unsharp outer boundary for the convective heat transport. It also contributes significantly to the energy transport in the interior (top left panel of Fig. 3), where convection never carries close to 100% of the luminosity. In the optically thick part the stratification is slightly far from radiative equilibrium and the entropy jump from the photosphere to the layers below is fairly large (entropy s in Figs. 3 and 4). The He ii/He iii ionization zone is visible in the entropy structure as a small minimum near the surface before the normal steep entropy decrease. The opacity peak at around T = 13   000 K just below the photosphere causes a very steep temperature jump (temperature T and opacity κ in Fig. 4), which is very prominent on top of upflow regions (opacity κ in Fig. 3). This causes a density inversion (density ρ in Figs. 4 and 8), which is a sufficient condition of convective instability resolved, the entropy drop occurs in a very thin layer, while the smearing from the averaging procedure over nonstationary up- and downflows leads to the large apparent extent of the mean superadiabatic layer.

The local radiative relaxation time scale in the photosphere is much shorter than a typical hydrodynamical time scale (time in Fig. 3). Numerically, the radiative energy exchange terms are quite stiff and prevent large deviations from a radiative equilibrium stratification. These terms enforce very short time steps (about 600 s per individual radiation transport step) and are responsible for the relatively high computational cost of this type of simulations. Local fluctuations in opacity, temperature (source function), and heat capacity pose high demands on the stability of the radiation transport module. This is true to a lesser degree for the hydrodynamics module due to shocks, high Mach numbers, and small scale heights. A side effect of the steep and significant temperature jump is the increase in pressure scale height from small photospheric values to values that are a considerable fraction of the radius in layers just below the photosphere (pressure scale heigh H_P in Fig. 3). The non-gray model has systematically shorter time scales (Fig. 3, bottom row), which causes the increased smoothing efficiency of temperature fluctuations.

Fig. 4 Logarithm of temperature (top left panel), density (top right), opacity (bottom left) and entropy (bottom right) of a slice through the center from a snapshot of the RSG simulation st35gm03n07 in Table 1.

4.2. Temperature and density structures

The temperature structure for the 3D simulations is displayed in Fig. 5 as a function of the optical depth at λ = 5000 Å for all simulations used in this work. We also computed classical 1D, LTE, plane-parallel, hydrostatic marcs model atmospheres (Gustafsson et al. 2008) with identical stellar parameters, input data, and chemical compositions as the 3D simulations (Table 2). It should be noted that the MARCS models do not have identical input data and numerical scheme as 3D simulations but they are a good benchmark for a quantitative comparison. The optical depth scale in Fig. 5 has been computed with the radiative transfer code Optim3D (see Paper I, and Sect. 5.2) for all rays parallel to the grid axes with a x  − y positions  ≤ (x_center,y_center) ± (G ∗ 0.25) (i.e., a small square at the center of the stellar disk), where x_center and y_center are the coordinate of the center of one face of the numerical box, G is the number of grid points from Table 1, and 0.25 was chosen to consider only the central rays of the box.

The models using a frequency-independent gray treatment show stronger temperature fluctuations compared to the non-gray case (top right panel in Fig. 5), as was already pointed out in Ludwig et al. (1994) for local hydrodynamical models, because the frequency-dependent radiative transfer causes an intensified heat exchange of a fluid element with its environment tending to reduce the temperature differences. Moreover, the temperature structure in the outer layers of the non-gray simulation tends to remain significantly cooler than the gray case (up to  ~200 hundreds K, Fig. 6), close to the radiative equilibrium. The thermal gradient is then steeper in the non-gray model and this is crucial for the formation of the spectral lines (see Sect. 5.2). This effect has already been pointed out in the literature and a more recent example is Vögler (2004), who pointed out that the non-gray treatment of opacities is mandatory to compare solar magneto-convective simulations with the Muram code to the observations.

It is also striking all gray simulations largely diverge from 1D MARCS models, while the case of the thermal structure of the non-gray model is more complicated. The outer layers agree very well with a cool 1D MARCS model at 3430 K (Fig. 5, top right panel). At lgτ₅₀₀₀ ~ 1 (i.e., where the continuum flux comes from) the mean 3D temperature is warmer than the 1D: a hot 1D MARCS model at 3700 K is then necessary to reach a better agreement but, however, this 1D profile diverges strongly for lgτ₅₀₀₀ < 1.

Fig. 5 Thermal structures of the simulations in Table 1 as a function of the optical depth at λ = 5000 Å. Darker areas indicate temperature values with higher probability. The solid light blue curve is the average temperature, the red dashed line is the 1D MARCS model profile with surface gravity g (Table 2) and the green dotted-dashed line is the 1D MARCS model with surface gravity geff (see text). In the top right panel the 3D mean thermal profile is compared to a cool MARCS model at 3430 K for the outer layers and a hot one with 3700 K (magenta dotted-triple dashed line) for the continuum-forming region (see text).

While turbulent pressure is naturally included in the RHD simulations, it is modeled in 1D models assuming a parameterisation as (6)where ρ is the density, v_t the characteristic velocity, and β is approximately 1 and depends on whether the motions occur more or less isotropically. This pressure is measuring the force produced by the kinetic movements of the gas, whether caused by convective or other turbulent gas motions. In 1D models, assuming spherical symmetry, the equation of hydrostatic equilibrium is solved for (7)where (being r the radius, M the mass of the star and G the Newton’s constant of gravity), P_g = ℜρT/μ_mol the gas pressure (ℜ being the gas constant, μ_mol the molecular weight, and T the temperature), and P_rad the radiative pressure. Gustafsson et al. (2008) and Gustafsson & Plez (1992) showed that assuming P_rad = 0 like in 3D simulations, they could mimic the turbulent pressure on the models by using models with those effects neglected with an adjusted gravity: (8)Eventually, Gustafsson et al. (2008) chose to set v_t = 0 for all 1D models in their grid, advising those who would have liked a different choice to use models with a different mass or gravity, according to the recipe given in Eq. (8).

After averaging all necessary thermodynamical quantities in 3D simulations averaged over spherical shells, we computed g_eff using Eq. (8) with μ_mol = 1.3 (appropriate value for the atmosphere of RSGs). Figure 7 shows the behavior of g/g_eff for 3D simulations of Table 1. We used the effective surface gravity at the the radius position (red vertical lines in the figure) to compute new 1D models (see Table 2). Figure 8 shows the comparison of the 3D density structures with 1D models. The models with the new effective gravity agree better with the 3D mean profiles than the models with v_t set to zero (i.e., g_eff = g). The effects on the temperature structures of 1D models with different surface gravity and equal effective temperature is negligible as displayed in all panels of Fig. 5.

We conclude that 1D simulation of RSG stars must have a turbulent velocity that is not equal to zero, which accordingly to Eq. (8) gives the correct value of surface gravity. The effect of the turbulent pressure is a lowering of the gravity that can be calibrated on 3D simulations in the same way as in to Fig. 7.

5. Non-gray versus gray simulations

We used the 3D pure-LTE radiative transfer code Optim3D described in Paper I to compute spectra and intensity maps from the gray simulation st35gm03n07 and the non-gray simulation st35gm03n13 (Table 1). The code takes into account the Doppler shifts caused by convective motions. The radiative transfer equation is solved monochromatically using extinction coefficients pre-tabulated as a function of temperature, density, and wavelength. The lookup tables were computed for the solar chemical compositions (Asplund et al. 2006) using the same extensive atomic and molecular opacity data as the latest generation of MARCS models (Gustafsson et al. 2008). We assumed a zero micro-turbulence because the velocity fields inherent in 3D models are expected to self-consistently and adequately account for non-thermal Doppler broadening of spectral lines.

5.1. From one- to three- dimensional models: determination of micro- and macroturbulence

Radiation hydrodynamics simulations provide a self-consistent ab-initio description of the non-thermal velocity field generated by convection, shock waves and overshoot that are only approximated in 1D models by empirical calibrations. Thus, the comparison between 1D and 3D spectra requires the determination of ad hoc parameters like microturbulence (ξ_mic) and macroturbulence (ξ_mac) velocities that must be determined for the 1D spectra.

For this purpose, we used a similar method as Steffen et al. (2009). Standard 1D radiative transfer codes usually apply ξ_mic and ξ_mac broadening on spectral lines isotropically, independently of depth, and identically for both atomic and molecular lines. Hence, to achieve a good representation of the conditions throughout the 3D atmosphere, we selected a set of 12 real Fe I lines in the 5000–7000 Å range (see Table 3) such that (i) their equivalent width on the 3D spectrum range between few mÅ to 150 mÅ; and (ii) with two excitation potentials, low (~1.0 eV) and high (~5.0 eV). Using the 1D MARCS models with the same stellar parameters (plus the effective surface gravity, g_eff) of the corresponding 3D simulations (see Table 2), we first derived the abundance from each line using Turbospectrum (Plez et al. 1993; Alvarez & Plez 1998; and further improvements by Plez) from the 3D equivalent width. For the non-gray model we used the 1D MARCS model corresponding to the lower temperature (i.e., 3430 K) because its thermal structure is more similar to the 3D mean temperature profile in the outer layers where spectral lines form (see top right panel of Fig. 5). Note that we computed the 3D spectra using Optim3D with negligible expected differences (less than 5%; Paper I) with respect to Turbospectrum. The microturbulence velocity was then derived by requiring no trend of the abundances against the equivalent widths. The error on the microturbulence velocity was estimated from the uncertainties on the null trend. Once the microturbulence was fixed, the macroturbulence velocity was determined by minimizing the difference between the 3D and the 1D line profiles, and this for three profile types: radial-tangential, Gaussian and exponential. The error on the macroturbulence velocity is calculated from the dispersion of the macroturbulence from line to line. The reduced χ² is also determined on the best-fitting 1D to 3D line profiles. Hence, the automatic nature of the determination ensures an objective determination.

Fig. 6 Mean profiles (top) and temperature difference (bottom) of the gray and non-gray simulations from Fig. 5.

Fig. 7 Ratio between the g and geff from Eq. (8) for simulations st35gm03n07 (solid line) and st35gm03n13 (dotted line) in top panel, and st36gm00n04 (solid line) and st36gm00n05 (dotted line) in the bottom panel. We used the same snapshots as in Fig. 5 and 8. The vertical red lines are the approximate positions of the radii from Table 1.

Table 4 presents the derived values of the microturbulence and macroturbulence respectively for the 3D gray and non-gray models. The dispersion of the microturbulent velocities is reasonably low, so that the depth-independent microturbulence is a fairly good approximation. Note also that the microturbulence velocity for the gray and non-gray model are quite close. Although the macroturbulence profile varies from one to another, the reduced χ² show that all of them give about the same value. Note that the micro and macroturbulence standard deviations from the average velocities are systematically higher in the non-gray model. We found that the non-gray spectral lines generally show a more complex profile that affect the line fitting with 1D models, which causes higher standard deviations. An example for the non-gray model is shown in Fig. 9: the spectral lines look shifted and broadened. The line asymmetries are caused by the inhomogeneous velocity field emerging from optical depths where lines form. The lower temperature in the outer layers of non-gray model (Fig. 6) could cause lower pressure and density scale heights and consequently a faster density drops and stronger shocks. These shocks would complicate the shape of the non-gray spectral lines. The macroturbulence parameter used for the convolution of 1D spectra reproduces only partly the complex profile with larger differences in line wings.

Fig. 8 Density structures of the simulations in Table 1. Darker areas indicate temperature values with higher probability. The colored curves have the same meaning as in Fig. 5).

On a wavelength scale of the size of a few spectral lines, the differences between 3D and 1D models are caused by the velocity field in 3D simulations that affect the spectral lines in terms of broadening and shift. A 1D model does not take into account velocity fields, and ad hoc parameters for micro- and macroturbulence are used to reproduce this effect. The microturbulence and macroturbulence calibration results are pretty similar for gray and non-gray simulations even if the 3D-1D corrections of the resulting iron abundances are smaller for the non-gray model. The microturbulence values obtained in this work are comparable, albeit a little bit lower, with Carr et al. (2000) who found 1.7 ≲ ξ_mic ≲ 3.5 km s^-1 with a different method based on observed CO lines. Our values also generally agree to what Steffen et al. (2009) found for their RHD Sun and Procyon simulations (0.8 ≲ ξ_mic ≲ 1.8).

Fig. 9 Example of CO first overtone lines computed for the 3D non-gray simulation and the corresponding 1D MARCS model of Table 2. The ξmic turbulent velocity used is 1.5 km s-1 and the ξmac is 6.4 km s-1 with a radial-tangential profile (see Table 3).

5.2. Spectral energy distribution

Fig. 10 Top panel: spectral synthesis of TiO band transition A3Φ − X3Δ(γ) (top row) and spectral energy distribution (bottom row) for the same snapshots of the gray, st35gm03n07, and non-gray, st35gm03n13, simulations of Fig. 5. 3D spectra are compared to the corresponding 1D MARCS models (Table 2) with ξmic and ξmac from Table 4 (radial-tangential profile).

The cooler outer layers encountered in the convection simulations are expected to have a significant impact on the temperature sensitive features. This is the case for, in particular, molecular lines (e.g., TiO) and strong low-excitation lines of neutral metals, the line-formation regions of which are shifted outwards because of the lower photospheric temperatures (Gustafsson et al. 2008; Plez et al. 1992). The importance of opacity from molecular electronic transition lines of TiO are crucial in the spectrum of M-type stars. Figure 10 (top row) shows the synthesis of a TiO band. The strength of the transition depends on the mean thermal gradient in the outer layers (τ₅₀₀₀ < 1), where TiO has a significant effect, which is more pronounced in the non-gray model with respect to the gray one Fig. 6. A shallow mean thermal gradient weakens the contrast between strong and weak lines and this is visible in the molecular band strength, which is much stronger for the non-gray model.

The top right panel of Fig. 10 also shows that the TiO band strength of non-gray model is more similar to the cool 1D MARCS model at 3430 K than the hot one at 3700 K. This reflects the fact that the 3D mean thermal structure in the outer layers is very similar to 1D-3430 K model (Fig. 5, top right panel).

The approximation of gray radiative transfer is justified only in the stellar interior and it is inaccurate in the optically thin layers. So far, this approximation was the best we could do because the RHD simulations have been constrained by execution time. However, with the advent of more powerful computers, the frequency-dependent treatment of radiative transfer is now possible and, on a wavelength scale of one or more molecular bands, the use of this method is big step forward for a quantitative analysis of observations.

Fig. 11 Top row: map of the linear intensity in the TiO band transition A3Φ − X3Δ(γ) of Fig. 10. The range is [0; 3.5 × 105] for the gray simulation st35gm03n07 of Table 1 and [0; 4.5 × 105] erg cm-2 s-1 Å-1 for the non-gray simulation st35gm03n13. Bottom row: visibility curves obtained from the maps above. The visibilities were computed for 36 different azimuth angles 5° apart. Bottom right panel: visibility fluctuations with respect to the average value for the gray (solid line) and non-gray (dashed line).

The shape of the spectral energy distribution (SED) reflects the mean thermal gradient of the simulations. The absolute flux plots in the bottom row of Fig. 10 display that at lower resolution two important conclusions can be retrieved: (i) the spectrum based on the gray model is very different from the 1D spectra; (ii) the non-gray model shows that the SED, compared to the 1D MARCS model with the hot temperature (3700 K), displays almost no distinction in the infrared region, and a weaker distinction in molecular bands in the visible region and in the near-ultraviolet region. Using the prescriptions by Bessell et al. (1998) for the filters BVRJHK, we computed the expected colors for the gray and non-gray models as well as the corresponding MARCS models (Table 5). The non-gray 3D model shows smaller differences to the 1D model than the gray one. The stronger radiative efficiency in the non-gray case forces the temperature stratification to be closer to the radiative equilibrium than in the gray case, where convection and waves have a larger impact on it.

Josselin et al. (2000) probed the usefulness of the index V − K as a temperature indicator for Galactic RSGs, but Levesque et al. (2006), fitting TiO band depths, showed that V − K and V − R provides systematic higher effective temperatures in Galactic and Magellan clouds RSGs (an average +170 K and +30 K, respectively) than the spectral fitting with 1D MARCS models. Levesque et al. (2006) concluded that the systematic offset was probably caused by the limitations of static 1D models. Using the radius definition of Eq. (5), we find that the effective temperature from the SED of 3D non-gray spectrum is T_eff = 3700 K. However, in the optical the spectrum looks more like a 3500 K 1D model (judging from the TiO bandheads of Fig. 10, top right panel). Assuming that the 3D non-gray model spectrum is close to real star spectra and using V − K, this leads to a T_eff that is higher by about 200 K than the T_eff from TiO bands. This goes in the right direction considering the results of Levesque et al. (2006).

5.3. Interferometry: visibility fluctuations

Chiavassa et al. (2010b, 2011) showed that in the visible spectral region, the gray approximation in our modeling may significantly affect the intensity contrast of the granulation. Figure 11 (top row) shows the resulting intensity maps computed in a TiO electronic transition. The resulting surface pattern is connected throughout the granulation to dynamical effects. The emerging intensity is related to layers where waves and shocks dominate together with the variation in opacity through the atmosphere and veiling by TiO molecules.

The intensity fluctuations are linked with the temperature inhomogeneities, which are weaker in the thermal structure of the non-gray simulation (see Sect. 4.2 and Fig. 5). We analyzed the impact on the interferometric observables using the method described in Paper I. We computed visibility curves for 36 different rotation angles with a step of 5° from the intensity maps in Fig. 11. In the plots, we introduced a theoretical spatial frequency scale expressed in units of inverse solar radii (). The conversion between visibilities expressed in the latter scale and in the more usual scale is given by . The resulting visibility curves are plotted in the bottom row of Fig. 11 together with the one σ visibility fluctuations, F, with respect to the average value (). The visibility fluctuations of the non-gray model are lower than the gray model ones for spatial frequencies lower than the first null point (approximately 0.0075 ), as a consequence the uncertainty on the radius determination (the first null point of the visibility curves is sensitive to the stellar radius, Hanbury Brown et al. 1974) is smaller:  ~10% in the gray model versus  ~4% in the non-gray one. At higher frequencies, the visibility fluctuations are higher in the non-gray model between 0.0075 and 0.002 , then the non-gray model fluctuations are systematically lower between  ~0.002 and 0.004 . However, it must be noted that after 0.03  (corresponding to 33 R_⊙, i.e.,  ~4 pixels), the numerical resolution limit is reached and visibility curves can be affected by numerical artifacts (Paper I and Sect. 6).

6. Study on the increase of the numerical resolution of simulations

We presented two simulations with approximately the same stellar parameters but with an increase of the numerical resolution. We analyzed the effect of increasing the resolution from 255³ grid points of simulation st36gm00n04 to 401³ grid points of st36gm00n05 (Table 1). Owing to the long computer time needed for the high-resolution model, the better way to proceed is to compute a gray model with the desired stellar parameters first, then the numerical resolution can be increased, and eventually the frequency-dependent non-gray treatment can be applied. With the actual computer power and the use of the OpenMP method of parallelization in CO⁵BOLD, this process may take several months.

Fig. 12 Maps of the square root intensity (to better show the structure in the higher resolution simulation) of st36gm00n04 and st36gm00n05 (Table 1) at 5000 Å together with a visibility curve extracted for a particular angle (solid curve) and the the same visibility’s angle after applying a [3 × 3] median smoothing (dashed curve). The intensity range in the maps is [0; 547.7] erg cm-2 s-1 Å-1.

We computed intensity maps at 5000 Å. This spectral region, as for the TiO band, is interspersed by a granulation pattern with high-intensity contrast. Figure 12 displays the intensity maps permeated by a variety of convection-related structures; while the large structures are visible in both images, the higher resolution model resolves the details much better and smaller-scale structures appear within the large ones (as already shown in a previous work, Chiavassa et al. 2006). Artifacts owing to Moire-like patterns show up in the maps of the higher numerical resolution simulation. They are composed of thin ripples up to 40% brighter than the surrounding points and have already been pointed out in Paper I. They are caused by high opacities, large variations in opacity, and significant changes in optical depth between two successive columns in the numerical box (see κ in Fig. 3). To overcome this problem, we computed visibility curves from these images and compared the visibility curves for one particular projected baseline from the raw image, and after applying a median [3 × 3] smoothing that effectively erases the artifacts. The visibility curves in Fig. 12 show that the signal from the low-resolution st36g00n04 simulation is stronger than st36g00n05 for spatial frequencies higher than about 0.07 . A weaker signal at a particular frequency means better spatial resolution: therefore very small structures (corresponding to roughly  ≤  15–20 R_⊙) are better resolved in st36g00n05.

Moreover, the visibility curves are mostly affected by the numerical artifacts for frequencies higher than  ~0.05 (corresponding to 25 R_⊙, i.e.,  ~5 pixels) for st36g00n04 and 0.09 (corresponding to 11 R_⊙, i.e.,  ~4 pixels) for st36g00n05. The increase of numerical resolution reduces only partly the problem in the intensity maps and therefore the cosmetic median filter must be still applied. This will not affect the visibilities at lower frequencies, which are the only ones observable in practice with modern interferometers.

Apart for the problems in the Moire-like patterns of intensity maps, the number of convection related surface structures seems to visually increase and possibly change size passing from 255³ to 401³ grid points (Fig. 12). In addition to this, the bottom right panel of Fig. 3 shows that going from 255³ to 401³ grid points the resolution of the temperature jump does not change so much. This means that the numerical resolution limit has not been reached yet. The principal limitation to the computation of an even more extended numerical simulation is the computer power. A more complete study concerning the impact of numerical resolution on the number and size of granules is necessary and will be addressed in a forthcoming paper, where a series of simulations with the same stellar parameters and increasing numerical resolution will be analyzed.

7. Conclusions

We presented the radiation hydrodynamical simulations of red supergiants with the code CO⁵BOLD together with the new generation of models based on a more sophisticated opacity treatment, and higher effective temperature and surface gravity.

The actual pattern and the mean temperature structure of RSGs simulations are influenced by (i) differences in the relative efficiency of convective and radiative energy transport and by the efficiency of radiation to smooth out temperature fluctuations; (ii) the optical depth of the upper boundary of the convection zone; (iii) the extent of convective overshoot.

The main difference between RSG and Sun-like granulation (except for geometrical scales) comes from the altered role of radiation: it is much more efficient in transporting energy in an RSG than in the Sun. This strongly influences the stratification and keeps it closer to the case of radiative equilibrium than its inefficient counterpart in the deeper layers of the Sun. Moreover, the strong entropy gradient causes a high buoyancy and renders the convective motions, which compete with the radiative energy transport, more violent in a relatively limited region over which the sub-photospheric entropy jump occurs (compared to a typical size of a convective element). Eventually, higher velocities are accompanied by higher pressure fluctuations and a stronger influence of shockwaves on the photosphere.

The most important improvement of our work is the relaxation of the assumption of gray opacities through multigroup frequency-dependent opacities. The non-gray simulation shows (i) a steeper mean thermal gradient in the outer layers that affect strongly the molecular and line strengths that are deeper in the non-gray case; (ii) that the general shape of the spectrum of the 3D non-gray simulation is similar to the 1D model, while the 3D gray simulation is very different. Hence, we concluded that the frequency-dependent treatment of radiative transfer with 5 bins represents a good approximation of the more complex profile of 1D based on  ~110 000 frequency points. Moreover, the lower temperature fluctuations of the non-gray model, caused by the intensified heat exchange of a fluid element with its environment, affect the surface intensity contrast and consequently the interferometric observables, which in turn reduces the uncertainty on the stellar radius determination.

We also showed that 1D models of RSGs must include a turbulent velocity calibrated on 3D simulations to obtain the effective surface gravity that mimics the effect of turbulent pressure on the stellar atmosphere.

We provided a calibration of the ad hoc micro- and macroturbulence parameters for 1D models using the 3D simulations, which have a self-consistent ab initio description of the non-thermal velocity field generated by convection, shock waves and overshoot. We found that the microturbulence velocity for the gray and non-gray model are quite close and the depth-independent microturbulence assumption in 1D models is a fairly good approximation, even if the 3D–1D corrections of the resulting iron abundances are smaller for the non-gray model. We also assessed that there is no clear distinction between the different macroturbulent profiles needed in 1D models to fit 3D synthetic lines; nevertheless, we noticed that micro and macroturbulence standard deviations on the average velocities are systematically larger in the non-gray model, which shows more complex line profiles than the gray simulation. This may be caused by stronger shocks in the outer layers of the non-gray models where the pressure and density scales heights are smaller.

While simulations with more sophisticated opacity treatment are essential for the subsequent analysis, models with a better numerical resolution are desirable to study the impact of the numerical resolution on the number and size of granules. We provided a first investigation showing that the number of convection related surface structures seems to increase and change size passing from 255³ to 401³ grid points, therefore the numerical resolution limit has not been reached yet. Future work will focus on a series of simulations with the same stellar parameters and increasing numerical resolution.

Acknowledgments

The ULB team is supported in part by an Action de recherche concertée (ARC) grant from the Direction générale de l’Enseignement non obligatoire et de la Recherche scientifique – Direction de la Recherche scientifique – Communauté française de Belgique. A.C. is supported by the F.R.S.-FNRS FRFC grant 2.4513.11, T.M. by the F.R.S.-FNRS FRFC grant 2.4533.09. A.C. thanks the Rechenzentrum Garching (RZG) and the CINES for providing the computational resources necessary for this work. A.C. also thanks Pieter Neyskens and Sophie Van Eck for enlightening discussions. B.F. acknowledges financial support from the Agence Nationale de la Recherche (ANR), the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, and the “École Normale Supérieure” (ENS) of Lyon, France, and the Istituto Nazionale di Astrofisica/Osservatorio Astronomico di Capodimonte (INAF/OAC) in Naples, Italy.

References

All Tables

All Figures

	Fig. 1 Gray intensity on one side of the computational cube from the initial sequence of the model st35gm00n05 in Table 1. The axes are in solar radii. The artifacts caused by the mismatch between the spherical object and the Cartesian grid become less evident with time passing.
In the text

	Fig. 2 Luminosity, temperature, and radius as a function of time for the simulations of Table 1: from left to right columns st35gm03n07, st35gm03n13, st36gm00n04, and st36gm00n05. The bottom panels are the ratio between turbulent pressure and gas pressure for different snapshots. The red vertical lines in all the panels is the approximative position of the radius from Table 1.
In the text

Fig. 3 Some spatially average quantities (i.e., spherical shells) for a snapshot of st35gm03n07 (Table 1). The red vertical line is the location of the radius. Top row, left panel: luminosities as function of radius in solar radii: continuous curve is the total luminosity; the dashed curve is the luminosity owing to pressure work; the dotted curve is the luminosity of kinetic energy; and the dash-dotted curve is the radiative luminosity. Top row, center: opacity. Top row, right: specific entropy. Central row, left panel: pressure scale height. Central row, center: characteristic radiative (continuous line) and hydrodynamical (dotted line) time scales. Central row, right: gravitational potential Φpot. Bottom row, left panel: random column opacity and radiative time scale for a given snapshot of the gray model st35gm03n07 (black) and non-gray model st35gm03n13 (red). The five colored curves correspond to the opacity groups. The radiative time scale is estimated for temperature fluctuations with spatial scales on the order of the grid size. Bottom row, center: characteristic radiative time scale for the gray (black) and non-gray (red) models. Bottom row, right: temperature profiles for the st36g00n04 (2553 grid points, black) and st36g00n05 (4013 grid points, red).

In the text

	Fig. 4 Logarithm of temperature (top left panel), density (top right), opacity (bottom left) and entropy (bottom right) of a slice through the center from a snapshot of the RSG simulation st35gm03n07 in Table 1.
In the text

Fig. 5 Thermal structures of the simulations in Table 1 as a function of the optical depth at λ = 5000 Å. Darker areas indicate temperature values with higher probability. The solid light blue curve is the average temperature, the red dashed line is the 1D MARCS model profile with surface gravity g (Table 2) and the green dotted-dashed line is the 1D MARCS model with surface gravity geff (see text). In the top right panel the 3D mean thermal profile is compared to a cool MARCS model at 3430 K for the outer layers and a hot one with 3700 K (magenta dotted-triple dashed line) for the continuum-forming region (see text).

In the text

	Fig. 6 Mean profiles (top) and temperature difference (bottom) of the gray and non-gray simulations from Fig. 5.
In the text

	Fig. 7 Ratio between the g and geff from Eq. (8) for simulations st35gm03n07 (solid line) and st35gm03n13 (dotted line) in top panel, and st36gm00n04 (solid line) and st36gm00n05 (dotted line) in the bottom panel. We used the same snapshots as in Fig. 5 and 8. The vertical red lines are the approximate positions of the radii from Table 1.
In the text

	Fig. 8 Density structures of the simulations in Table 1. Darker areas indicate temperature values with higher probability. The colored curves have the same meaning as in Fig. 5).
In the text

	Fig. 9 Example of CO first overtone lines computed for the 3D non-gray simulation and the corresponding 1D MARCS model of Table 2. The ξmic turbulent velocity used is 1.5 km s-1 and the ξmac is 6.4 km s-1 with a radial-tangential profile (see Table 3).
In the text

	Fig. 10 Top panel: spectral synthesis of TiO band transition A3Φ − X3Δ(γ) (top row) and spectral energy distribution (bottom row) for the same snapshots of the gray, st35gm03n07, and non-gray, st35gm03n13, simulations of Fig. 5. 3D spectra are compared to the corresponding 1D MARCS models (Table 2) with ξmic and ξmac from Table 4 (radial-tangential profile).
In the text

Fig. 11 Top row: map of the linear intensity in the TiO band transition A3Φ − X3Δ(γ) of Fig. 10. The range is [0; 3.5 × 105] for the gray simulation st35gm03n07 of Table 1 and [0; 4.5 × 105] erg cm-2 s-1 Å-1 for the non-gray simulation st35gm03n13. Bottom row: visibility curves obtained from the maps above. The visibilities were computed for 36 different azimuth angles 5° apart. Bottom right panel: visibility fluctuations with respect to the average value for the gray (solid line) and non-gray (dashed line).

In the text

	Fig. 12 Maps of the square root intensity (to better show the structure in the higher resolution simulation) of st36gm00n04 and st36gm00n05 (Table 1) at 5000 Å together with a visibility curve extracted for a particular angle (solid curve) and the the same visibility’s angle after applying a [3 × 3] median smoothing (dashed curve). The intensity range in the maps is [0; 547.7] erg cm-2 s-1 Å-1.
In the text