1 Introduction

The structure and dynamics of the solar photosphere is crucially determined by the mass and energy transport processes taking place across the solar surface. For the study of the solar convection and its phenomenological manifestation on the visible surface of the Sun, the solar granulation, numerical simulations not only complement observational data but also serve as a means of their own by providing complete and almost continuous information in 3-D of the physical state and the dynamics which otherwise often has to be drawn indirectly from observations. The physics of the layers surrounding the solar surface is rather involved; below the surface the opacity is sufficiently large so that the local adiabatic gradient is exceeded by the temperature gradient needed for radiative-diffusive energy transport, turning the fluid convectively unstable. This process is primarily described by mixing length theory (e.g. Cox and Giuli 1968; Spruit 1974) that proved to be successful at determining the average energy transport (Cattaneo et al. 1991), while neglecting e.g. the dynamical modifications to the hydrostatic equilibrium near the surface (MacGregor 1991).

At the solar surface and above, along with a rapid decrease of the opacity, radiation becomes the primary energy transport mechanism, although a continuing mass flow overshooting into convectively stable regions still characterizes the lower photosphere. Ionization and molecular dissociation processes taking place in the near surface layers considerably affect the internal energy and the equation of state. The extent to which numerical models reflect reality has improved significantly with increasing numerical accuracy and level of physical detail since the 1970s, one such enhancement being the consideration of non-gray radiative transport.

The photosphere and the subjacent thin \(4~\mbox{Mm}\) wide layer of the upper convection zone that are covered in the simulation studied here constitute a region of almost negligible radial extent, accounting for no more than 0.64% of the solar radius. Yet it is these layers that feature remarkably rich dynamical processes driven by steep vertical gradients of temperature, pressure, and density, the latter falling off by 5 and 4 orders of magnitude, respectively, see Fig. 1.

Fig. 1
figure 1

Left: Vertical profile of horizontally averaged pressure, density, and temperature for a given model time step as a function of geometric depth \(x\) as obtained from the ANTARES model based on the ATLAS 9 package (Kurucz 1970). Error-bars indicate the quantities’ variation at a given geometric depth. Right: Comparison to data from the energy-balance model atmosphere No. C of the quiet Sun of Fontenla et al. (1993). The steep temperature gradient at the surface that is visible in both model atmospheres is due to the large temperature sensitivity of the dominant \(\mathrm{H}^{-}\) opacity (e.g. Stein and Nordlund 1998)

Full size image

We present the state-of-the-art radiation hydrodynamics (RHD) code ANTARES (A Numerical Tool for Astrophysical RESearch), applied to the study of the near surface convection and the photosphere of the Sun. In spite of the broad range of applicability, ranging from the modeling of photospheric turbulence (Muthsam et al. 2007) to Cepheid pulsation (Muthsam et al. 2011), the code has not received much attention in the solar physics community so far. While, for the time being, the code is restricted to the modeling of the quiet Sun, an MHD upgrade is foreseen in the near future that will allow us to study the modifications to the dynamics and to the energy transport due to the photospheric magnetic field and to extend its applicability to magnetoactive regions at the solar surface. Further recent developments of the code are summarized in Blies et al. (2015) and include the consideration of two-component flows (Zaussinger 2010), a parallel multigrid solver for the 2-D non-linear Helmholtz equation (Happenhofer et al. 2013), and a generalization to solve the Navier-Stokes equations on curvilinear grids (Grimm-Strele et al. 2014). While the over many years of development fully matured ANTARES code is characterized by its elaborate numerical schemes and its stability, it is by far not the only simulation project of similar scope of applicability. Of the many other codes we can only exemplarily mention a few such as the 3-D MHD code of Nordlund and Galsgaard (1997) that lead the way in many respects and by which it was first possible to model the solar atmosphere from the photosphere up to the corona (Gudiksen and Nordlund 2002). Of late the code CO5BOLD, developed at the Kiepenhauer Institut (e.g. Freytag et al. 2010) has gained a lot of attention in connection with the observation of rapidly rotating magnetic structures, so-called magnetic tornadoes, in chromospheric simulation data that serve as an energy channel into the corona (Wedemeyer-Böhm et al. 2012). It not so long ago obtained an optional generalization from RHD to MHD (Schaffenberger et al. 2005) and, in conjunction with the related spectral analysis package Linfor3D (Gallagher et al. 2016) is also capable of computing the emergent spectrum, allowing for a more direct comparison to observations. The interaction of magnetic fields with radiative convection and further magnetic activity in the photosphere, such as the formation and dynamics of magnetic flux tubes are intensely studied with the MURaM code (e.g. Vögler et al. 2005), that is a collaborative project of the Max-Planck Institute for Solar System Research and the University of Chicago. Finally, the young computational heliophysics 2-fluid code JOANNA that is developed at the Maria Curie-Skłodowska University in Lublin, Poland already revealed promising results concerning the simulation of spicules, showing that triggering pulses steepen into upward propagating shocks whereas the chromospheric cold and dense plasma lags behind such shocks during its rise into the corona with a mean speed of \(20\ldots25~\mbox{km}/\mbox{s}\) (Kuźma et al. 2017). This is neither a ranking list nor should it be considered to be anything but incomplete, its sole purpose here being to give an impression of the industriousness of this particular research area and of the vast efforts that are put into the examination of the complex energy transport mechanisms in the solar atmosphere with i.a. the intention to finally fully understand the heating of the solar corona.

The high resolution of the ANTARES model photospheres allows further detailed studies of smaller-scale dynamical phenomena such as quiet-Sun jets that have been reported from the IMaX instrument on board of the SUNRISE stratospheric balloon telescope (Martínez Pillet et al. 2011; Borrero et al. 2010), jet-like vortex tubes that have been observed with the New Solar Telescope (NST) by Yurchyshyn et al. (2011), or fine-scaled MHD phenomena occurring in MHD simulations (Kitiashvili 2014) that cause intense dynamic interactions between the surface and the chromosphere and which may be responsible for acoustic wave excitation and quasi-periodic flow eruptions, see e.g. Kitiashvili et al. (2013) and references listed therein. The present RHD model has already revealed rotating plasma jets (Lemmerer et al. 2016), which seem to be triggered by turbulent convection. Our current investigations (Lemmerer et al. 2016) suggest that horizontal flows around rotating jets may trigger MHD kink waves (e.g. Zaqarashvili et al. 2007) or torsional Alfvén waves (Fedun et al. 2011; Shelyag et al. 2013; Zaqarashvili et al. 2013) that may propagate in magnetic flux tubes that are anchored in the photosphere and thereby transport photospheric energy into the chromosphere. First approximations showed that the observed kinetic energy flux may excite waves carrying an energy flux of \(\sim 7 \times10^{7}~\mbox{erg}\,\mbox{cm}^{-2}\,\mbox{s}^{-1}\), supposing only a 1% energy transfer to the waves. This wave energy flux is about one order higher than necessary to compensate for energy losses from the quiet-Sun chromosphere and two orders higher than needed to heat the quiet corona (Withbroe and Noyes 1977). Already a 10% wave energy dissipation into heat would be sufficient to heat the chromosphere and corona. Further investigations based on soon available ANTARES RMHD (radiation magnetohydrodynamics) photospheric models will be crucial to test these first assumptions and to determine modifications due to the involvement of the photospheric magnetic field.

Correlation analysis of characteristic parameters, which is the focus of the present preliminary study, is an approved method for examining the solar granulation dynamics and photospheric stratification and dates back to as early as the 1950s, when e.g. Stuart and Rush (1954) found a correlation between velocity and brightness—a result that pointed to the no longer surprising fact that bright granules are associated with convective matter upflow. Durrant and Nesis (1982) found that this correlation extends from the continuum to a height of 300 km beyond that layer, corresponding to a “pure convective component” of the velocity field scaling about \(3''\) in size, while the velocity field on smaller scales appeared to be more turbulent. Height dependent correlations were analyzed by Gadun et al. (2000) and references listed therein to study the vertical photospheric structure. The former found the top convection zone to reach 20 to 50 km into the photosphere, from where an onset of convective overshoot in stable regions into a height of 150 to 170 km was observed. Beyond that zone the columnar photospheric structure is still maintained due to the influence of convective pressure variations up to a height of about 250 to 300 km, a zone referred to as transition layer and being characterized by its pronounced inversion of temperature fluctuations. A breakdown of the columnar structure takes place at a height of 300 km, from where oscillations govern the dynamics of the photospheric medium. These findings based on coherence analysis confirmed prior results from spectral observations (Nesis et al. 1988; Karpinsky 1990).

The remainder of this paper is organized as follows: In Sect. 2.1 the fundamental equations of radiation hydrodynamics are introduced. We further outline important numerical methods and discuss their advantages for the reliability of the model results. Also boundary and initial conditions as applied in the particular model presented here are addressed. In Sect. 2.2 some basic techniques for the analysis of the simulation data are presented. For the study of the photospheric structure we first examine the height variation of typical model quantities’ horizontal distributions in Sect. 3.1 and discuss some of their correlations directly at the solar surface. Relative fluctuations of the temperature, gas pressure, and density give important insights into the diverse dynamics of the different photospheric layers and are discussed in Sect. 3.2. A fuller picture of the photospheric stratification is gained from studying the correlation of characteristic model quantities locally as well as by a two-point correlation, where one quantity is fixed at the surface while the other one is varied with height, which is the objective of Sect. 3.3. Finally, in Sect. 4 we discuss the model results and present possibilities for further investigation.

2 Methods

2.1 Model description

The study at hand is based on a simulation of the photospheric matter by the code ANTARES (Muthsam et al. 2007, 2010). For the quiet Sun the dynamics of the solar surface layers can be described by the equations of radiation hydrodynamics, i.e. the continuity equation

$$ \frac{\partial\varrho}{\partial t} + \nabla\cdot(\varrho\boldsymbol{u}) = 0, $$
(1)

Euler’s equation of momentum balance

(2)

and an energy balance

$$ \frac{\partial e}{\partial t} + \nabla\cdot\bigl(\boldsymbol{u} (e + P) - \boldsymbol{u} \cdot\boldsymbol{\tau} \bigr) = \varrho(\boldsymbol{g} \cdot \boldsymbol{u}) + Q_{\mathrm{rad}}, $$
(3)

complemented by the equation of state closing the conservation laws (1) to (3).Footnote 1 Here \(\varrho\) denotes the mass density, \(\boldsymbol{u} = u \boldsymbol{e}_{x} + v \boldsymbol{e}_{y} + w \boldsymbol{e}_{z}\) is the flow velocity, is the stress tensor with the kinetic gas pressure and the viscous stress tensor , where \(\mu\) is the dynamic molecular viscosity. Furthermore \(\boldsymbol{f} = \varrho\boldsymbol{g}\) denotes an external force density with \(\boldsymbol{g}\), the local gravitational acceleration, is the momentum current density tensor, is the total kinetic energy density made up of the convective kinetic- and the internal energy density, and \(Q_{\mathrm{rad}}\) represents radiative sources.

In the scope of RHD finally the frequency-dependent and time-independentFootnote 2 radiative transfer has to be considered. The radiation transfer equation

$$ \hat{\boldsymbol{r}} \cdot\nabla I_{\nu}= \varrho\kappa_{\nu}(S_{\nu} - I_{\nu}) $$
(4)

with spectral intensity \(I_{\nu}\), source function \(S_{\nu}\), and material opacity \(\kappa_{\nu}\) is solved for the upper \(\sim1~\mbox{Mm}\) of the computational domain. It is linked to the energy balance equation via the radiative heating rate \(Q_{\mathrm{rad}}= - \int(\nabla\cdot \boldsymbol{F}_{\nu}) \mathrm{d} \nu\), accounting for the energy exchange between the gas and the radiation field, where \(\boldsymbol{F}_{\nu}\), the frequency \(\nu\)-dependent radiative energy flux is the spectral intensity \(I_{\nu}\) integrated over the solid angle \(\varOmega\) into that energy is radiated along unit vector \(\hat {\boldsymbol{r}}\), i.e. \(\boldsymbol{F}_{\nu}= \int I_{\nu}(\hat{\boldsymbol{r}}) \hat{\boldsymbol{r}} \mathrm{d} \varOmega\). By solving the equations of radiation-hydrodynamics, the code simulates convection and full (i.e. non-gray) radiative transfer in local thermal equilibrium (LTE).Footnote 3 The full radiative flux \(\boldsymbol {F}_{\mathrm{rad}}\) is found by integrating \(\boldsymbol{F}_{\nu}\) over frequency using non-gray opacities with \(N=4\) bins, applying the frequency binning method developed originally by Nordlund (1982) and used extensively thereafter. At a sufficient depth, typically \(\tau \gtrsim100\) (Steiner et al. 1997), the heating rate can be computed according to the diffusion approximation \(Q_{\mathrm{rad}} = \nabla\cdot(\kappa\nabla T)\).

The temporal integration of the model equations uses second- or third-order Runge-Kutta methods with weighted essentially non-oscillatory (WENO) finite volume schemes. The code is based on a high-resolution finite-volume method that can treat turbulence by adopting local mesh refinement. Essentially, finite volume schemes are based on interpolation of discrete data using polynomials; fixed stencil interpolations work well for sufficiently smooth problems but introduce oscillations near discontinuities, whose amplitudes do not decay in the course of mesh refinement. Whereas traditional remedies such as the introduction of an artificial viscosity or the application of limiters to discard such oscillations have obvious drawbacks, ENO (originally developed by Harten and Osher 1987; Harten et al. 1987) and Weighted ENO schemesFootnote 4 (Liu et al. 1994; Jiang and Shu 1996) are based on a nonlinear adaptive procedure to determine the locally smoothest stencil, whereby the crossing of discontinuities in the interpolation procedure can be avoided by and large (Shu 1998). Further information concerning the implementation of WENO in ANTARES can be found in Kupka et al. (2012), Mundprecht et al. (2013), Zaussinger and Spruit (2013), and Grimm-Strele et al. (2014).

In horizontal directions, periodic boundary conditions for all quantities are applied. The current model has however been improved lately w.r.t. the applied vertical boundary conditions. In principle, boundary conditions should influence the flow kinematics at the boundaries as little as possible. By originally inhibiting any vertical convective matter and energy flux through the lower and upper boundaries by setting the vertical flow velocity component to zero, \(u|_{\mathrm{top}} = u|_{\mathrm{bot}} = 0\), mass and energy conservation had been ensured. While such closed boundary conditions have the advantage of simplicity and high stability, they disturb the velocity field in an undesirable way (Robinson et al. 2003) and are in general prone to non-physical reflections of waves and shocks and had been justified as the optical depth unity layer had only been weakly influenced (Mundprecht et al. 2013). Also the energy fluxes have been found to be influenced by forcing zero fluid velocity to an extent of about two pressure scale heights above the lower end of the computational domain (Grimm-Strele et al. 2015). The current model remedies these deficiencies by applying open boundary conditions at the top and at the bottom, thereby allowing for convective mass and energy in- and outflows. This is of special significance at the lower boundary where most of the energy transport is due to advection and kinetic energy flux. The lack of convective matter flux into the bottom layer in prior models had to be compensated for by an appropriate artificial radiative source term. The boundaries are made transmissive for waves by constantly extrapolating velocities according to

$$ \frac{\partial u}{\partial x} = \frac{\partial v}{\partial x} = \frac {\partial w}{\partial x} = 0, $$
(5)

thereby increasing the long-term stability of the model (Grimm-Strele et al. 2015). The latter study gives a detailed description of the recent numerical implementations of boundary conditions in the ANTARES code and also compares the effects of closed and several open boundary conditions.

For initial conditions a solar atmosphere between layers of \(4350~\mbox{K}\) at the top and \(20\,000~\mbox{K}\) at the bottom, respectively is considered. The horizontal momentum densities \(\varrho v\) and \(\varrho w\) are slightly disturbed before the system is evolved in time.

The simulation time exceeds 5 hours with a time step of \(\Delta t = 8.7~\mbox{s}\) between snapshots. The model domain is covered by a Cartesian lattice \(\varOmega= \{(i,j,k)|i,j,k \in\mathbb {N}_{0}; 0 \leq i \leq404; 0 \leq j,k \leq510\}\) ranging from the top of the photosphere to \(\sim4~\mbox{Mm}\) into the upper convection zone and covering a field of \(18~\mbox{Mm} \times18~\mbox{Mm}\) in horizontal direction. The considerably high spatial horizontal resolution of \(\Delta y = \Delta z = 35.3~\mbox{km}\) exceeds the resolution of recent observational data, e.g. images of the 1 m Swedish Solar Telescope (SST) that, neglecting the loss of resolution capacity due to seeing conditions, has a theoretical diffraction limit of e.g. \(\lambda/D \approx 0.14 ''\) in \(\mbox{H}\alpha\) (Antolin and Rouppe van der Voort 2012), corresponding to \(\sim 100~\mbox{km}\) on the solar disc, but also of the 1.6 m Solar Telescope at Big Bear Observatory with a diffraction limit of \(\approx0.07''\) at 500 nm (Goode et al. 2002) and of the 1.5 m GREGOR telescope at the Observatorio del Teide, Tenerife which resolves spatial scales on the solar surface down to \(\approx 50~\mbox{km}\) (Puschmann et al. 2012). The vertical dimension of the model grid is even better resolved with \(\Delta x \approx 11.0~\mbox{km}\), enabling us to accurately examine the photospheric stratification. The top of the grid is associated with zero optical depth. With increasing opacity, the solar surface is found to be located on average \(\approx 600~\mbox{km}\) below.

2.2 Data analysis techniques

A first insight into the structure of the solar photosphere is gained from analyzing relative fluctuations \(\delta Q = (Q - \langle Q \rangle)/\langle Q \rangle\) of a quantity \(Q\), such as the thermodynamic state variables temperature \(T\), gas pressure \(P\) and mass density \(\varrho \). Here \(\langle Q \rangle\) denotes a horizontal average of the particular quantity in question, evaluated separately in upflow- \(U = \{(i,j,k)|u_{i jk} \leq0\}\) and downflow regions \(D = \{(i,j,k)|u_{i jk} > 0\}\) with vertical flow velocity \(u = \boldsymbol{u} \cdot\boldsymbol{e}_{x}\). The vertical profiles of these fluctuations are studied after having been finally averaged horizontally over up- and downflow regions and over time.

The linear dependency between two matrix-valued quantities \(A_{jk}\) and \(B_{jk}\) evaluated at horizontal surfaces for a given vertical level \(i\) is described by the linear correlation coefficient

$$ \rho_{i} = \frac{\sum_{j,k} (A_{i jk} - \bar{A}_{i} ) (B_{l jk} - \bar{B}_{l} )}{\sqrt{ (\sum_{j,k} ( A_{i jk} - \bar{A}_{i} )^{2} ) ( \sum_{j,k} (B_{l jk} - \bar{B}_{l} )^{2} )}}, $$
(6)

which for \(l = i\) measures the linear dependence at the same depth, also referred to as local one-point correlation. In general, \(l \neq i\), which is also referred to as two-point correlation and from which changes in the columnar photospheric structure can be analyzed (Gadun et al. 2000). Note that e.g. \(\bar{A}_{i}\) simply denotes the horizontal mean of \(A_{i jk}\) for fixed \(i\). Both kinds of methods are used extensively in the following to measure the respective local and non-local correlations between fluctuations of temperature \(\delta T\), pressure \(\delta P\), mass density \(\delta\varrho\), opacity \(\delta\kappa\) and vertical fluid flow component \(\delta u\).

The solar surface corresponding to \(\ln\tau= 0\), i.e. optical depth unity, serves as a reference level for the two-point correlations defined above. The boundary at the upper vertical level of the grid \(x_{i} = 0\) is associated with a zero optical depth \(\tau_{0,jk} = 0\ \forall j,k\), from where by a post-processing numerical integration top-down a columnar estimate for the optical depth is obtained,

$$ \tau_{i jk} \approx\tau_{i-1,jk} + \int_{x_{i-1}}^{x_{i}} \kappa_{jk}(x) \varrho_{jk}(x) \, \mathrm{d} x \quad (i \geq1). $$
(7)

Figure 2 shows the logarithmic run of \(\langle \tau\rangle\) with geometric depth, indicating a maximum optical depth of \(\sim10^{8}\) at the lower end of the lattice inside the convection zone corresponding to a depth of \(\approx3.9~\mbox{Mm}\) below the visible surface. From this estimate, Eq. (7), the mean level of \(\tau= 1\) as obtained by averaging over columns and time is found to be located at \(\langle x|_{\tau= 1} \rangle\approx 602.8~\mbox{km}\). The evaluated \(\tau\)-unity isosurface and its horizontal average for a snapshot in a \(2 \times2~\mbox{Mm}^{2}\) subfield are shown in Fig. 3. As the \(\mbox{H}^{-}\) opacity increases with temperature, the corrugated isosurface juts out above the averaged geometric depth level in the hot upflow regions, while in the intergranular lanes deeper layers are observed. The rms-variation of optical depth unity is \(\approx 34.8~\mbox{km}\), a value comparable with the optical depth corrugation found from the solar granulation models of Stein and Nordlund (1998). Since at the constant geometric depth level temperatures in the hot upflow regions are measured deeper down in even hotter layers and higher up in the cooler downflow regions as compared to the visible surface, a much broader temperature- and intensity range is to be expected. This is illustrated in the left panel of Fig. 4, where the temperature profile along a horizontal and a \(\tau\)-unity slice placed across two minor and a major granule of \(\sim1~\mbox{Mm}\) size clearly shows a significantly larger temperature (and accordingly also intensity-) variation when evaluated at the averaged constant geometric depth level rather than at the \(\tau\)-unity isosurface itself.

Fig. 2
figure 2

Profile of the temporally and spatially averaged optical depth \(\tau(x)\) as a function of geometric depth as obtained by numerical integration of Eq. (7)

Full size image
Fig. 3
figure 3

Temporal snapshot of \(x|_{\tau=1}(y,z)\) (red) and the horizontal average \(\langle x|_{\tau=1} \rangle_{y,z}\) (black). Higher layers, i.e. smaller values of \(x\) correspond to the granular brightness field, while impressed regions are associated with the intergranular lanes, where the solar surface is observed at deeper layers

Full size image
Fig. 4
figure 4

Left: Temperature and relative intensity runs at a horizontal slice at constant geometric depth corresponding to the horizontal mean of optical depth unity (solid) and along a slice lying at the \(\tau=1\) isosurface (dashed). Right: Temperature distribution \(T(x,y)\) (K) and \(\tau=1\) isosurface at a vertical cross section

Full size image

3 Results

3.1 Stratification of the photosphere

Figure 5 shows the height variation of typical model quantities’ horizontal distributions observed at a given snapshot in the photosphere. As we can see from the temperature and intensity images in the first and second row, respectively, the temperature maxima at the surface and below are located above the granules, while the gas in the intergranular lanes is some thousand degrees cooler. Higher up in the photosphere—the chosen level of 275 km in the rightmost column of the figure corresponds to a maximum temperature difference in up- and downflows in the middle photosphere as will be discussed in the following sections—the temperatures are inverted due to buoyancy breaking such that the gas above granules is significantly cooler than the coalescing downdrafts in the intergranular lanes. In contrast, the observed intensity in the middle photosphere does not mirror its image at the surface: granules are brighter throughout the photosphere. The intensity distributions peak at the borders of the granular upflows and in particular in regions bordering merging intergranular lanes. As discussed above, temperatures and intensities vary significantly stronger at a given geometric depth than at a surface of constant optical depth. Comparing the intensity and velocity distributions (third row) shows that the highest intensity regions within granules correspond also to the maximum upward directed flow velocities. Downdrafts are observed to be highest where at least two intergranular lanes meet and the cooled flows of several granules converge. Again, the range of values is here slightly higher at surface level-constant geometric depth, but maximum values for downflow velocities are found at constant optical depth unity, where downflow regions are located in considerably deeper layers. The density and pressure distributions are shown in rows 4 and 5. Below the surface coalescing downflows are significantly denser than the granular upflows and are associated with a lower pressure. While pressure and density images are mirrored in subphotospheric layers, relative density- \(\delta\varrho\) and pressure fluctuations \(\delta P\) become of comparable size above the transition layer once radiative equilibrium is established. Here, the energy exchange of the perturbations with the surroundings can be considered isothermal. Above the transition layer we observe a gradual breakdown of the columnar structure. The mean and extreme values of temperature, intensity, gas density, intensity and vertical flow velocity are summed up in Table 1, each for granular and intergranular regions on constant geometric depth and on the \(\tau\)-unity isosurface, respectively.

Fig. 5
figure 5

Distributions of temperature, intensity, vertical flow velocity, gas density, and gas pressure evaluated at characteristic height levels in the photosphere. The variation of the horizontal temperature distribution \(T(y,z)\) at the visible surface is much smaller than the one at the zero-height level as hot granules are observed higher up due the temperature-increasing \(\mbox{H}^{-}\) opacity while cooler gas in the intergranular lanes is observed at considerably deeper layers

Full size image
Table 1 Mean and extreme values of various quantities at the solar surface for a given time step. Each quantity is given separately for granular and intergranular regions on a horizontal plane corresponding to horizontally averaged optical depth unity and on the optical depth unity isosurface evaluated from Eq. (7)
Full size table

In the following we will put this qualitative description of the spatial variation of model quantities with the granular brightness field as well as with height on a more quantitative basis in terms of a correlation analysis. We start our analysis of the model photosphere by reproducing the most obvious correlation between the vertical flow velocity and the brightness at the solar surface that was found in observations already in the mid-twentieth century (e.g. Stuart and Rush 1954). From the leftmost panel of Fig. 6 it is apparent that bright granular areas at the \(\tau\)-unity isosurface are associated with a negative, i.e. upward vertical fluid flow velocity and vice versa for the darker gas in the intergranular lanes. The correlations of the vertical velocity with temperature (middle panel) and gas density (right panel) show that at the surface the upflows are significantly hotter but less dense compared to the intergranular downdrafts. These trivial dependencies already reveal that the correlation of the vertical flow with the density is weaker than its correlation with the brightness, just as that the allocation of low density gas with upflows is a good one but not unflawed. This incident was already apparent from column 3 of Fig. 5. The wide dispersion shows that low gas densities are to be found in downflows as well and the strongest downdrafts with \(u \approx8~\mbox{km}/\mbox{s}\) even exclusively exhibit low gas densities in contrast with the overall trend.

Fig. 6
figure 6

Correlation of intensity, temperature, and gas density with the radial flow velocity at the visible surface \(\tau=1\)

Full size image

3.2 Relative fluctuations of thermodynamic variables in two-component representation

Relative fluctuations \(\delta Q\) of some quantity \(Q\) were evaluated according to the procedure described in Sect. 2.2. The convectively unstable subphotosheric layers are characterized by divergent temperature fluctuations in up- and downflows that are most prominent only 90 km and 70 km, respectively below the photosphere, Fig. 7. As one moves further upwards into the photosphere, temperature fluctuations in up- and downflows reverse their sign at a height of \(\approx 130~\mbox{km}\). This inversion that peaks at \(\approx275~\mbox{km}\) can be attributed to the buoyancy breaking effect explaining the inversion via two mechanisms: First, as the upflowing matter overshoots into the convectively stable region it is bound to lose energy via radiation due to a strong decrease in opacity. Second, as downflows are coalescing into the intergranular lanes, they become compressed and heated. Going still further upwards, the temperature of upflows increases again as these upper layers become compressed by oscillations. These height levels were found to vary depending on the employed photospheric model: The first reversal of temperature for instance is located in a height range from 100 km based on the mechanical-radiative energy balance model of the solar granulation by Musman and Nelson (1976) up to 170 km as observed from the 2-D RHD model of Gadun et al. (1999). A slight variation is found also in relative height levels such as the distance between the two temperature reversal points that in the latter 2-D model falls short of the distance observed from our modeling by 20%.

Fig. 7
figure 7

Relative temperature fluctuations averaged over up- (red) and downflows (black) as well as time as a function of height from subphotospheric layers to the convectively stable adjacent photosphere. The inset is a zoom of the photospheric region showing the subdivision into layers of convective overshoot above the surface, a transition layer, and oscillating layers in the upper photosphere

Full size image

Figure 8 shows that absolute pressure fluctuations are significantly greater than absolute density fluctuations \(|\delta P| > |\delta\varrho|\) in the lower photosphere and become of roughly equal size above the first reversal of temperature fluctuations \(|\delta P| \approx|\delta \varrho|\), indicating radiative equilibrium above that level, while directly at the reversal points their concurrence is exact, \(|\delta P|=|\delta\varrho|\).

Fig. 8
figure 8

Relative fluctuations of pressure (solid) and density (dashed) averaged over upflows (red), downflows (black) and time

Full size image

Temporally and horizontally averaged vertical velocities evaluated separately for granular and intergranular flows are shown in Fig. 9. Error-bars indicate the temporal variability of the horizontal mean of the quantity in consideration. The vertical velocity field does not invert across the photosphere. Maxima of the vertical flow velocity in up- and downdrafts as listed in Table 1 are averaged out and hence are not reflected in this graph. Downdrafts on average have higher velocities throughout the photosphere as well as in subphotospheric layers. Both, up- and downdrafts reach their maximum flow velocity just below the photosphere, where also the temperature separation of the two flows is most pronounced. In absolute values the average flow velocities of up- and downdrafts decrease from the surface across the lower photosphere toward the top of the transition zone and rise again in the higher oscillatory layer. As was pointed out already by Schwarzschild (1948), acoustic waves above granules proceed essentially without dissipation into the higher photosphere and the perturbations do not affect the temperature there, such that the energy flux in acoustic waves \(\varrho \boldsymbol{u}^{2} c_{\mathrm{s}}\) with sound speed \(c_{\mathrm{s}}\) can be considered constant, explaining the rise of the average flow velocity with decreasing density in the oscillatory layers. This requires flow velocities well below a critical value \(u^{2} < u_{\mathrm{crit}}^{2} = \beta P/\varrho\) with some fraction \(\beta\) of the order of 0.1. While by and large the same argument also applies to the downflows, one should not leave unmentioned that locally fast downdrafts exceed the critical flow speed and even shocks are found sparsely scattered in the intergranular lanes.

Fig. 9
figure 9

Absolute values of the vertical velocity in up- and downflows for the full computational domain

Full size image

3.3 Local and two-point model quantity correlations

We first discuss local or one-point linear correlations between various model quantities before turning to (relative) two-point correlations for which one quantity is fixed at the solar surface while the second one is varied with height.

The one-point correlation between the vertical velocity- and temperature fluctuations \(\rho(\delta u, \delta T)\) pictured in the upper panel of Fig. 10 again reflects the temperature reversals discussed before. As expected, a strong correlation is found in the subphotospheric convective layers and an anticorrelation in the most part of the photosphere due to the overcooling of the photospheric matter in optically thin layers. We also observe from this figure that the anticorrelation of the vertical velocity- and gas density fluctuations \(\rho(\delta u,\delta\varrho)\), characteristic for convectively unstable regions, reverses its sign at a height of \(\approx35~\mbox{km}\) up from where a pursuing positive correlation is found. The gas density and pressure are strongly correlated throughout the photosphere, see the lower panel of Fig. 10.

Fig. 10
figure 10

One-point correlations (averaged over time) as a function of height in the photosphere. Error-bars indicate the temporal variation of these correlation functions. Top: Correlation of fluctuations of the vertical velocity with temperature (red) and gas density (black). Bottom: Correlation of fluctuations of the gas density with gas pressure (red) and the vertical velocity with gas pressure

Full size image

A high correlation of temperature fluctuations at the solar surface with temperature fluctuations at varying heights, \(\rho(\delta T, \delta T_{i})\) is found from subphotospheric layers to the lower photosphere as shown in the upper panel of Fig. 11. As expected, the correlation steeply drops from the solar surface (where it is of course exact) upwards until it becomes negative below the temperature reversal of up- and downflows. Throughout the region of temperature reversal the anticorrelation is strongest and is found to decrease again above the second temperature reversal point yet without turning positive again. In subphotospheric layers the hotter and brighter gas of the granulation cells is less dense as can also be seen by comparing the temperature and density images at \(h=-100~\mbox{km}\) in Fig. 5, why not surprisingly a strong anticorrelation \(\rho(\delta T, \delta\rho_{i}) \approx-0.6\) is found where the temperature separation of up- and downflows has reached its maximum level few ten kilometers below the surface, see lower panel of Fig. 11. As the temperature reverses from here, the correlation function becomes positive very fast, reaching a high level of \(\rho(\delta T, \delta\rho_{i}) \approx0.6\) in the transition region between thermal convection and the oscillating layers. Here, hotter and brighter gas is found in the intergranular lanes, cf. also the temperature image at \(h = 275~\mbox{km}\) in Fig. 5 which almost mirrors the temperature distribution at \(h = -100~\mbox{km}\). As the temperature separation in up- and downflows diminishes above that layer again, so also decreases the correlation in the above regions where the columnar structure is no longer present.

Fig. 11
figure 11

Top: Two-point correlations of temperature fluctuations (fixed at \(x|_{\tau=1}\)) and running temperature and pressure fluctuations, respectively. Bottom: Same for temperature fluctuations and running density as well as vertical velocity fluctuations

Full size image

The gas pressure is higher in granular upflows throughout subphotospheric and photospheric regions resulting in an entirely positive correlation \(\rho(\delta T, \delta P_{i}) \geq 0.2\) (black curve in the leftmost panel of Fig. 11). Of course, throughout and above the transition layer, where due to a radiative equilibrium \(\delta P \approx\delta\rho\), the correlation functions of temperature fluctuations with \(\delta P_{i}\) and \(\delta\varrho_{i}\), respectively show a similar run with height.

From the images of the vertical velocity in Fig. 5 it is apparent, that \(u\) is always negative above granules (corresponding to upflows) and positive in the intergranular lanes. In the oscillatory layers however that clear allocation becomes increasingly blurred as the columnar structure breaks down. Quantitatively this is shown by the correlation function \(\rho(\delta T, \delta u_{i})\) (black curve in the lower panel of Fig. 11). The high correlation in subphotospheric layers and the lower photosphere drops notably to insignificant values above the transition layer.

Finally, we discuss the two-point correlation of the intensity fluctuations at the surface with the mean opacity fluctuations at varying heights, \(\rho(\delta I, \delta\kappa _{i})\), cf. the red curve in Fig. 12. It is striking that up to the surface this correlation function almost coincides with the correlation of intensity fluctuations with temperature fluctuations, which, as was already argued by Gadun et al. (2000), is due to the ionization of hydrogen which is strongly temperature dependent as \(\mbox{H}^{-}\) ions primarily cause the absorption of photons here. In higher layers where due to the still high temperatures basically ionized metals are mainly accountable for the absorption of photons, the fluctuations of the opacity are no longer sensitive to temperature fluctuations. The authors’ claim based on their 2-D RHD model that up from here \(\rho(\delta I, \delta\kappa _{i})\) is closely following the correlation function \(\rho(\delta I, \delta P_{i})\) could not be reproduced with our recent model, Fig. 12.

Fig. 12
figure 12

Correlation of intensity fluctuations (fixed at \(x|_{\tau =1}\)) with running opacity-, pressure- and temperature fluctuations

Full size image

Altogether these data provide the necessary information for outlining the overall vertical structure of the photosphere (see also Table 2): Due to the rapid decrease in opacity, see right panel of Fig. 2, radiative cooling quickly gains in importance as one proceeds from subphotospheric layers across the \(\tau=1\) isosurface. Following the argument of Gadun et al. (2000), we can interpret the height level where \(\rho(\delta u,\delta\varrho)\) turns positive as the top of the thermal convection zone which thus extends to some ten kilometers above the solar surface. While the adjacent layer still exhibits positive and negative temperature fluctuations in up- and downflows, the further increasing positive correlation \(\rho(\delta u, \delta\varrho)\) indicates that this layer is no longer convectively unstable although thermal convective upflows overshoot into this region up to a height of \(\approx 130~\mbox{km}\), where a radiative equilibrium is established and temperature fluctuations reverse their sign. The columnar structure of thermal convection persists up to a height of \(\approx275~\mbox{km}\) where the temperature inversion is most pronounced. The adjacent upper photosphere is governed by acoustic oscillations and the columnar structure is no longer observable. Comparing these findings with the ones of pioneering RHD simulations of the solar granulation developed by Musman and Nelson (1976) or the 2-D model of Gadun et al. (2000), we find a good agreement of the qualitative runs of correlation functions; finally we were able to update the structural division levels based on our up-to-date numeric model.

Table 2 Photospheric stratification classified by its dynamically distinguished layers
Full size table

4 Discussion

We introduce the radiation hydrodynamics code ANTARES that we applied to the study of the solar granulation and that has not yet received much attention in the Solar Physics community. We used correlation analysis to examine the vertical stratification of the photosphere and determined height levels subdividing the photosphere in layers that exhibit characteristic dynamics of their own: The subphotospheric layers up to a height of \(\approx35~\mbox{km}\) above the solar surface were found to be convectively unstable. Convective upflows overshoot further into the lower photosphere into a height of \(\approx 130~\mbox{km}\), where temperature fluctuations in up- and downflows coincide and become exactly zero. The overlying layer is a transition region between the convective and oscillatory regimes. Within its roughly 145 km extension the horizontal distribution of the gas temperature mirrors the one at the photospheric footpoints and below. That inversion peaks at the top of this layer. Further up the columnar structure of the photosphere gradually breaks down as one proceeds into the oscillation-controlled higher photosphere. This rough schematic structure confirms findings from previous models (e.g. Gadun et al. 2000; Musman and Nelson 1976) and spectral observations, (e.g. Nesis et al. 1988; Karpinsky 1990), while with the present model some further accuracy to the sensitively model-dependent height levels is added.

The WENO scheme implemented in ANTARES avoiding oscillatory solutions at discontinuities which otherwise occur due to the interpolation of discrete data in finite volume methods is particularly useful for the ongoing study of shocks which are observed in the intergranulum of our model photospheres and for the study of acoustic oscillations in the scope of RHD. We intend to further investigate photospheric wave excitation and propagation by the application of wavelet and wavepacket analysis and to quantify the associated energy transfer through the photosphere.

As this RHD-code is currently heavily under development with an imminent RMHD upgrade to be released, we intend to soon present further model results and focus on photospheric, small-scale, intergranular rotating plasma jets that have been detected and studied in our RHD simulations (Lemmerer et al. 2016) but whose supposed contribution to the chromospheric and coronal heating via the generation of MHD kink waves or torsional Alfvén waves relies on testing our assumptions by studying equivalent RMHD model photospheres.