NON-GREY RADIATIVE TRANSFER IN THE PHOTOSPHERIC CONVECTION : VALIDITY OF THE EDDINGTON APPROXIMATION

Abstract

The aim of this study is to describe the physical processes taking place in the solar photosphere. Based on 3D hydrodynamic simulations including a detailed radiation transfer scheme, we investigate thermodynamic structures and radiation fields in solar surface convection. As a starting model, the initial stratification in the outer envelope calculated using the solar calibrations in the context of the standard stellar theory. When the numerical fluid becomes thermally relaxed, the thermodynamic structure of the steady-state turbulent flow was explicitly collected. Particularly, a non-grey radiative transfer incorporating the opacity distribution function was considered in our calculations. In addition, we evaluate the classical approximations that are usually adopted in the onedimensional stellar structure models. We numerically reconfirm that radiation fields are well represented by the asymptotic characteristics of the Eddington approximation (the diffusion limit and the streaming limit). However, this classical approximation underestimates radiation energy in the shallow layers near the surface, which implies that a reliable treatment of the non-grey line opacities is crucial for the accurate description of the photospheric convection phenomenon.

1. INTRODUCTION

The structure of low-mass stars is closely connected to convective energy transport in their outer envelopes. Especially, the superadiabatic layers (SAL) at the top of the convection zone (CZ) is a highly turbulent region where the energy transport radically changes from convection to radiation. In order to accurately describe the convective motions in this radiating fluid, a reliable treatment of both convection and radiation is required. Although the underlying relations of a radiating fluid have been well-established, the non-linear characteristics of the governing equations is extremely unwieldy. For this reason, the standard stellar model depends on numerous classical approximations such as the mixing length theory, the macro-/micro- turbulence parameter, the Eddington T − τ relation and the grey atmospheres. Through empirical tuning of the physical parameters among observational data, 1-dimensional (1D) stellar models have successfully provided the global properties of stars. However, convection and radiation itself is a 3-dimensional (3D) phenomenon, and therefore, the 1D stellar model can not properly describe the photospheric convection.

With advances in computing power, numerical simulations of the 3D radiation-hydrodynamics (RHD) have enabled a direct computation of the surface convection using the fewest number of conventional approximations (Nordlund 1982; Chan & Wolff 1982; Stein & Nordlund 1989). However, a detailed treatment of radiative transfer in the hydrodynamic simulations is still computationally demanding, and so few multi-dimensional models for the Sun and solar-type stars have been investigated (Nordlund& Dravins 1990; Steffen et al. 1989; Kim et al. 1995; Freytag et al. 1996; Robinson et al. 2005). Recently, considering a wide range of stellar evolutionary states and fundamental parameters, a fine grid of the RHD models have been compiled (Bach & Kim 2012; Trampedach et al. 2014a,b; Tanner et al. 2014). Moreover, the recent solar chemical compositions have been determined by spectral line synthesis based on the 3D RHD simulations (Asplund et al. 2005, 2009). The significance of the RHD simulation is to provide a reliability test of the classical approximations that have been inevitably prescribed in the 1D stellar models. Firstly, the mixing length parameter (Böhm-Vitense 1958) that presents a convective efficiency is highly uncertain. On the basis of multi-dimensional simulations, a numerical calibration of this mixing length parameter has been suggested (Chan & Sofia 1987; Kim et al. 1996; Ludwig et al. 1999; Trampedach & Stein 2011). In general, the mixing length parameter substantially affects the thermodynamic structure in the upper convection zone of low-mass stars, which implies that atmospheric boundary should be defined at the same time. Through intensive model computations, an absolute calibration of the mixing length parameter has been provided (Trampedach et al. 2014b), and then it was combined with a more realistic T − τ relation (Trampedach et al. 2014a). If the classical approximation is properly calibrated with a large variety of stellar parameters, a more reliable 1D model can be achieved with cheaper computational effort.

Interactions between radiation and matter are the key processes in the photospheric convection of stars. These interactions are statistically characterized by a macroscopic coefficient, that is, opacities of the medium. Due to the frequency dependency of radiation fields, the grey approximation using mean opacities oversimplifies the radiative transfer problem especially in the shallow region of stars. Moreover, isotropic thermal diffusion combined with the grey atmospheres can not account for the non-local characteristics of radiation fields in the optically thin layers. In their 2D solar simulations, Vögler et al. (2004) examined a nongrey radiative transfer by comparing the opacity distribution function with the multigroup opacities. They suggested that at least 3-5 binned opacities are required to imitate the total radiative heating computed with the opacity distribution function. In this study, we investigate the 3D hydrodynamic surface convection including a detailed radiative calculation. Adopting opacity distribution functions, a fully nongrey treatment has been applied to the 3D hydrodynamical simulations. The ultimate aim of this study is to (1) reproduce the surface convection on the basis of the 3D radiationhydrodynamics, (2) improve radiative transfer computations considering a detailed non-grey atmospheres, (3) determine the validity of the classical Eddington approximation.

 

2. HYDRODYNAMICS

Hydrodynamic simulation is one of the main tools in the studies on the astrophysical gaseous systems. The solar granule is an extremely fast-moving region driven by the photospheric convection. To realistically describe the solar surface convection, the following requirements should be satisfied :

Assuming no eletromagnetic interactions and no chemical reactions, the governing equations for the fully compressible Newtonian fluid can be written as

in a conservative form. Here, E = e + ρv2/2 is the total energy density containing the internal energy density e = ρcVT. Thermodynamic quantities ρ, v, P and g denote density, velocity, pressure and gravitational acceleration. The hydrodynamic terms are coupledwith the radiative heat energy Qrad for the forward cooling and the backward heating process. In cases of low-mass stars, momentum exchange due to radiation pressure can be safely neglected. The diffusive flux f denotes the heat diffusion defined by f = −k∇T, which involves k is the thermal conductivity of the gas. Generally, the momentum and energy transfer causes an internal friction, which is captured by the viscous stress tensor

where the strain rate tensor σij = (vi,j + vj,i)/2, the Kronecker delta δij, the dynamical viscosity coefficient μ, the dilatational viscosity coefficient λ, and the bulk viscosity ζ (= λ+2μ/3). Under conditions of high temperature and pressure inside stars, the fluid is regarded as a perfectmonatomic gas. Then, the bulk viscosity becomes zero, and the Stokes relation 2μ + 3λ = 0 is satisfied.

Large-Eddy simulation (LES) has become one of a most promising methodology for reproducing turbulent flows. Preferentially, it postulates that most energy is transported by resolved large eddies that define the geometry of the flow, while small eddies are implicitly accounted for by a subgrid-scale (SGS) parameterization (Smagorinsky 1963; Deardorff 1970). Then, the eddy viscosity of the resolved scale is defined as

where the colon denotes contraction of the strain rate tensor, the sub-grid scale eddy coefficient cμ, and the local mesh size Δ = (ΔxΔy)1/2Δz. When an unphysical shock takes place, the eddy viscosity is momentarily controlled by 1 + C · (∇ · v)2 maintaining numerical stability.

 

3. RADIATIVE TRANSFER

3.1. Non-Grey Atmospheres

The super adiabatic layer at the top of the outer envelope is a transition region of the energy transport. For a realistic description of radiant flow of the photospheric convection, a detailed treatment of radiative transfer is required. In this study, we investigate a non-grey radiative transfer in the 3D hydrodynamical medium. Compared with our previous work (Bach & Kim 2012), two main enhancements have been applied to the transfer computations. Firstly, a set of transfer equation is directly solved using the operator splitting method. Secondly, adopting the opacity distribution functions (Kurucz 1995), we consider a non-grey atmospheres instead of the mean opacities.

In a plane parallel atmospheres, the radiative transfer equation is written by

where Iμν is the specific intensity, the source function Sν, the directional cosine μ, and the optical depth τν. A general transfer problem can be established by a set of non-linear, integro-differential equations with respect to spatial, angular, and wavelength of radiation fields. From the angular integrals, three moment equations are given by

where Jν is the mean intensity, Hν the Eddington flux, Kν the K–integral. Meanwhile, the frequency-integrated mean intensity is defined as

where φν is the normalized absorption profile. Assuming complete redistribution and isotropic scattering, the thermal source function S is defined as

where B is the Planck function and ϵ is the photon destruction probability per scattering. Then, the transfer problem reduces to two equations of Equation (4) and (9) (Mihalas 1978).

By introducing a linear operator

a formal solution for the transfer problem is

where I is the identity matrix. The linear operator Λ is called the ordinary Λ-operator that relates radiation fields to the local source function. Due to the angular-/frequencydependence of the radiation fields, however, this complete linearization requires a high computational cost. Particularly, it converges slowly in computations of the stellar atmospheres (Mihalas 1978; Hubeny 2003).

Instead, for a suitably chosen approximate operator Λ∗,

the iteration can be written by

and the acceleration term (I − αΛ∗)−1 becomes diagonalized (Cannon 1973a,b). This is the well-known accelerated lambda iteration (ALI) method (Hubeny 2003). Indeed, this numerical acceleration scheme greatly improved the rate of convergence in our transfer computations. The iteration was set to terminate when the maximum change in the source function reaches |ΔS/Sn| ≲ 10−5. Because the grid spacing is relatively sparse and only few points exist near the surface, the vertical grid mesh has been re-discretized using a logarithmic scale according to the optical depth. For deep layers (τ ≥ 104), the asymptotic relation of thermal diffusion has been incorporated.

In order to consider the effects of the line transfer, we employed the opacity distribution function provided by Kurucz (1995). The basic idea of the opacity distribution function is that absorption coefficients may be represented by a continuous probability distribution functions that do not vary radically over the frequency intervals. In practice, the number and the width of the frequency intervals are empirically determined. Particularly, we adopt recently updated opacity distribution functions (Sbordone et al. 2004) and subroutines of the ATLAS9 stellar atmospheric code (Kurucz 1995, 1996; Castelli et al. 1997) in our transfer computation. Line opacities were treated with scattering coefficients (Sν ≈ Jν) and continuous opacities were computed separately. By summing over all opacity groups, total absorption coefficients were determined.

3.2. The generalized Eddington Approximation

In a static medium, the radiative term can be separated from the heat equation. If a small perturbation is imposed on an isothermal medium, the radiative relaxation of thermal energy can be solved exactly (Spiegel 1957). The diffusion approximation implies that the radiative heat can be treated by thermal diffusion. In deep layers, absorption coefficients are assumed to be independent of frequency (κν ≡ κ) in a completely redistributed medium. For isotropic specific intensity, the complete solution of the transfer equation is called the Eddington approximation. Although the Eddington approximation oversimplifies the transfer problem, it has provided valuable understanding of energy balance in the stellar atmospheres. Considering an anisotropic source function, Unno & Spiegel (1966) generalized the Eddington approximation,

Based on this generalized transfer equation, the surface convection phenomenon of the Sun (Kim et al. 1996; Robinson et al. 2003), Procyon A (Robinson et al. 2004), and subgiant stars (Robinson et al. 2005) has been studied. More recently, Bach & Kim (2012) investigated dynamical properties of turbulent motions generated from three most recent solar mixtures (Grevesse & Sauval 1998; Asplund et al. 2005, 2009). The main goal of this study is to improve the radiation schemes in our 3D hydrodynamic code and to evaluate the classical Eddington approximation.

 

4. SIMULATIONS

4.1. Initial Configurations

The starting model for numerical simulations were obtained from the standard solar calibration adopting Grevesse & Sauval (1998) chemical compositions, which are in good agreement with the helioseimological inversion (Basu & Antia 2008). The initial stratification of the thermodynamical quantities such as internal energy, density, and temperature were calculated using the Yale stellar evolutionary code (YREC). In our solar calibrations, the mixing length parameters and helium abundance are adjusted such that the standard sun has the radius R⊙ and the luminosity L⊙ of the present-day value.

Microscopic chemical diffusion is a slow migration of small particles due to gravitational settling and radiative levitation. Feeding heavy elements into the deep region, the microscopic chemical diffusion can modify internal constituents of stars. Especially, in cases of low-mass stars, even a small rate of diffusion can considerably affect local physical properties (equation of state and opacities) during their long main-sequence lifetime. Adopting the standard diffusion coefficients (Bahcall & Loeb 1990; Thoul et al. 1994), we include the microscopic diffusion of helium (Y) and heavy elements (Z) in the evolutionary computations. Taking the solar age tage ∼ 4.55 Gyr, the surface chemical composition is determined to be (Xs, Ys, Zs) = (0.7407, 0.2442, 0.0169), while the initial abundance is estimated to be (X0, Y0, Z0) = (0.7090, 0.272, 0.0188) at the pre-main sequence (PMS) phase.

In the standard evolutionary model, the efficiency of the convective energy transport is usually quantified by the classical mixing-length parameter αMLT = l/Hp (Böhm-Vitense 1958). Even a tiny variation of this parameter can significantly affect thermodynamic structures of the upper part of convection zone. In our solar calibrations, the solar mixing length parameter was determined to be α⊙ = 1.881, but it is unnecessary after starting the 3D simulations. Our evolutionary computation includes two sets of the Rosseland mean opacities : OPAL opacities (Iglesias & Rogers 1996, updated in 2006) in deep layers, and the low temperature opacities (Ferguson et al. 2005) near the surface. In the overlapping region (4.0 ≤ log Teff ≤ 4.5), these two opacities are averaged. We also adopt the OPAL equation of state (Rogers et al. 1996, updated in 2006). Stellar parameters and physical inputs in the solar calibration are summarized in Table 1. Thermodynamic variables and opacities were suitably converted into the initial configurations for the 3D simulations using dimensionless relations of the scaling factors (Chan & Wolff 1982). For the efficient growth of turbulence, a dynamical perturbation is imposed on the initial stratification.

Table 1The standard solar parameters

4.2. Numerical Simulations

The geometry of the computational domain is assumed to be a plane-parallel, closed box. The closed boundary condition guarantees conservation of mass, momentum and energy. However, this may lead to an unphysical shock and a boundary effect. When an unphysical shock takes place, it is controlled by the sub-grid eddy viscosity (Smagorinsky 1963). To minimize boundary effects, the top boundary is placed in the stable regions above temperature minimum layers (Tmin ∼ 4100 K) where almost all energy is transported by radiation. At the bottom boundary, a constant heat flux Fb = 6.3064×1010 erg sec−1cm−2 is supplied to the computing domain. If the flux entering fluctuates, it will be smoothed out in the strongly stratified layers over the Kelvin-Helmholtz timescale τKH. The numerical fluid is isolated with an impenetrable (stress-free) rigid top and bottom boundary. We also assume a periodic side-boundaries. The detailed numerical schemes of our large-eddy simulations are described in Appendix A. The size of the computing domain has physical dimension of 40002 × 3000 km that is divided into 1172 × 170 Eulerian cells. The time step Δt is adjusted to be 0.02 sec in the implicit phases and is set to be 0.01 sec in the fully explicit simulations. In general, the maximum size of time interval Δt is constrained by the local flow speed, the sound speed and the grid spacing. For stable convergence, our simulations have to satisfy a necessary condition known as the Courant-Friedrichs-Lewy (CFL) condition,

where the grid spacing Δ, the flow speed v, and the adiabatic sound speed Additionally, molecular and thermal diffusion among grids are treated by Δt ≤ Δ2/ν and Δt ≤ Δ2/κ. The height of computing domain corresponds to 10-12 pressure scale-height that completely contains the super adiabatic layers. The width of the box covers a number of typical solar granules. Our grid resolution is sufficiently fine to resolve turbulent convection near the solar surface (Bach & Kim 2012). A horizontal snapshot of the 3D structures is shown in Figure 1.

Figure 1.The horizontal slice of a 3D snapshot near the surface. The box size is 4,000 × 4,000 km with the resolution of 1172 ×170 using the Eulerian grid schemes. The computational domain of simulations covers several typical solar granules horizontally.

 

5. THE VALIDITY TEST OF THE CLASSICAL APPROXIMATIONS

For 200min of the real-time convection that sufficiently covers the solar convective turn-over timescale (τTO ∼ 45 min), the 3D structures have been computed by the fully explicit methods. In Figure 2, radiation fields in the vicinity of the solar surface are illustrated. The frequency-integratedmean intensities (Equation 8) have been computed at the height of Zh = 150, 0, -300, -1200 km from the solar surface. The horizontally averaged temperature at each height corresponds to T ∼ 4580, 5770, 12000, and 15000 K, respectively. Here, the frequency-integrated mean intensity is scaled by J0 = 1.0 × 1010 erg sec−1 cm−2. In order to verify dynamical properties, the vertical velocity distribution in the upper convection zone is presented in Figure 3. The root mean square (rms) velocities have been temporally averaged over the whole computational domain. In general, the solar granules are characterized by fast (narrow) downflow and slow (wide) upflow, which can be measured by the photospheric observations of the C-shape line bisector (Dravins 1987). Recent power spectra analysis of temperature gradient and velocity fields revealed that the fluctuations have vertical dimensions at granules (λ = 0.5-2.0 Mm), mesogranules (5-12 Mm), and supergranules (20-30 Mm) (Baran & Stodilka 2014; Hathaway et al. 2015). From the intensity distribution and the generated velocity contrast, our simulation successfully reproduces the typical solar granulations.

Figure 2.The frequency-integrated mean intensities are presented at the height of Zh = 150, 0, -300, -1200 km (from top to bottom panel) from the solar surface. The horizontally averaged temperature of each slice corresponds to T ∼ 4580, 5770, 12000, and 15000 K, respectively. The width of simulation domain is extended to 40002 km. is scaled by J0 = 1.0 × 1010 erg sec−1 cm−2.

Figure 3.The vertical distribution of the root mean square velocities in the upper convection zone is presented. Slow upflow in the wide regions (dotted) and fast downflow (solid) in the narrow regions characterize the typical solar granules. This velocity contrast implies that our numerical simulations successfully reconstruct the solar surface convection. The top axis denotes the vertical distance (not scaled) from the surface.

Adopting the fewest assumptions, a direct computation of the photospheric convection has become feasible (Chan & Sofia 1987; Stein & Nordlund 1989; Freytag et al. 1996). More specifically, numerical simulations can be used to test the oversimplified classical approximations. In previous studies (Chan & Sofia 1987; Kim et al. 1996; Ludwig et al. 1999), a numerical calibration of the classical mixing length parameter (Böhm-Vitense 1958) was suggested. Recently, the validity of the mixing length parameter was evaluated using fine model grids of the RHD simulation (Trampedach et al. 2014a,b). In addition, a more realistic T − τ relation defining the stellar boundary has been combined with the calibrated mixing length parameter (Trampedach et al. 2014b). The effect of metallicity on the surface convection zone was also investigated (Bach & Kim 2012; Tanner et al. 2014). However, hydrodynamic simulations including radiative transfer are still computationally expensive. Accordingly, the numerical verification of the classical approximations will provide an efficient way to treat a complex physical process inside stars.

In Figure 4, the vertical distribution of radiation fields is presented. The frequency-integratedmean intensities are averaged horizontally and temporally during the whole computational domain of the explicit simulations. Here, the frequency-integrated mean intensity is scaled by J0 = 1.0 × 1011 erg sec−1 cm−2. The dotted line denotes the generalized Eddington approximation (Equation 14), and the solid line is the accelerated Λ-iteration schemes using the opacity distribution functions. At both ends, the two transfer schemes are in good agreement. In the optically thick region (τ ≫ 1), the radiative energy is well represented by isotropic thermal diffusion defined as

where a is the radiation constant, c is the speed of light. This implies that the corresponding medium is completely redistributed and absorption coefficients are well approximated by the Rosseland mean opacities. In the optically thin layers (τ → 0), most radiant energy is released out of the surface, and only small amount of energy is absorbed in the atmospheres. Fromthe asymptotic properties, the source function distribution is well described by the two limits of the Eddington approximation : the diffusion limit and the streaming limit. However, in the shallow intermediate layers, the energy transport drastically varies from convection to radiation. The absorption coefficients decrease steeply, and the line opacities become important. Then, the Rosseland mean opacity as a harmonic mean considerably underestimates the strong line opacities, and the resultant radiation fields are overestimated (< 18%). This implies that the non-grey effect is crucial to understand the photospheric convection. However, a complete inclusion of the frequency dependent opacities is still unknown and computationally impractical. If the frequency spectrum is divided into a binned opacity group according to line strength, the absorbed energy can be represented by several numbers for each binned mean opacities (Nordlund 1982). As discussed in their 2D simulations, Vögler et al. (2004) suggested the use of multigroup opacities in multi-dimensional transfer problems. In their 2D solar simulations, Vögler et al. (2004) compared the multigroup opacities with the opacity distribution function. They suggested that at least 4-5 binned opacities are required to reliably account for radiative heating in the atmospheres.

Figure 4.The comparison test of the Eddington approximation and the lambda iteration in the hydrodynamic medium. The vertical distribution of the frequency-integrated mean intensities is presented. The dotted line represents radiation fields computed by the grey schemes of the Eddington approximation. The solid line denotes the result of the accelerated Λ-iteration method incorporated by the non-grey atmospheres of the opacity distribution function. Squares and crosses are grid points for the radiative transfer computations. The frequency-integrated mean intensity is scaled by J0 = 1.0 × 1011 erg sec−1 cm−2. The vertical lines denote typical optical depths τ = 2/3, 1.0, and 10 that are computed from the conventional T − τ relation.

 

6. SUMMARY

In order to describe radiative processes near the solar surface, we performed hydrodynamical simulations including the detailed transfer calculation. For initial configurations of the numerical simulations, the standard solar model (SSM) has been calibrated within the framework of the 1D stellar evolution theory. Using the initial stratifications of the thermodynamic variables, the 3D hydrodynamic simulations have been performed. When the numerical fluid thermally relaxed and the steady-state turbulence statistically grew, the 3D snapshots were compiled by fully explicit numerical schemes. Particularly, the transfer problem was solved directly using the accelerated Λ-iteration method based on the non-grey atmospheres of the opacity distribution functions.

From the size distribution of the typical granules and the velocity contrast of the numerical flow, the solar photospheric convection was successfully reproduced. Especially, the validity of the Eddington approximation was verified by comparing the vertical distribution of the radiation energy. From the asymptotic properties of the Eddington approximation, it well approximates the source function distribution in the optically thin and thick layers. Since the Rosseland mean opacity considerably underestimates strong lines, grey atmospheres overestimate the radiation energy in the SAL. In order to qualitatively estimate the effect of the detailed line-transfer, a non-grey treatment of opacities should be employed. To obtain a more realistic atmospheric structure, we are computing a finer grid simulations. In addition, a set of RHD models for low mass stars is in progress. A detailed description of radiation and convection will provide a deeper insight of physical processes inside stellar atmospheres.