© T. S. Bhatia et al. 2022

1 Introduction

The interpretation of data from stellar observations requires comparison against stellar models. Traditionally, these models have been 1D global models (Carbon & Gingerich 1969) that use formulations of mixing-length theory (MLT; Böhm-Vitense 1958). Later models have accounted for line-blanketing effects and used an opacity distribution function (ODF) approach to calculate opacities (Strom & Kurucz 1966). Among them, the MARCS code (Gustafsson et al. 1975, 2008), the ATLAS code (Kurucz 1979; Castelli & Kurucz 2003), and the PHOENIX code (Allard & Hauschildt 1995; Hauschildt et al. 1999) have enabled the calculation of synthetic stellar spectra with a detailed accounting for the relevant physics. These models have enabled, for example, the accurate determination of abundances and stellar evolution tracks, along with constraining the chemical evolution of galaxies (Edvardsson et al. 1993).

Nonetheless, convection is a 3D process and a phenomenological 1D approach is insufficient for characterizing the properties of granulation and plasma turbulence. Radiative aspects of granulation, in particular, are not fully captured by an MLT-like approach. The 1D models also require free parameters, such as the mixing length parameter, α_MLT, and turbulent velocities for line broadening calculations. Hence, 3D stellar atmosphere models become important for a realistic interpretation of stellar characteristics from observations. The first 3D simulations of solar convection were pioneered by Nordlund (1982); Stein & Nordlund (1998); Nordlund & Stein (1990). These simulations are realistic in the sense that they are directly comparable to solar observations: they reproduce granulation pattern and associated spectral line widths and asymmetries quite well (Asplund et al. 2000).

Early stellar models (Nordlund & Dravins 1990) showed the limitations of MLT-based models in accurately reproducing the near-surface temperature gradient, which affects radiative properties in the lower photospheres. Subsequently, various models have been constructed over a grid of effective temperature, surface gravity, and metallicity. Currently, the most comprehensive 3D grids (Magic et al. 2013; Ludwig et al. 2009) cover a wide range of stellar type on the HR-diagram, but they are purely hydrodynamic.

However, stellar convection is not a purely hydrodynamic process. Most cool stars are expected to have magnetic fields. Hence, a complete description of their photospheres should also take into account the effects of such fields. The best studied star in this context is the Sun. There is a rich variety of solar magnetic field-related phenomena ranging from sunspots and active regions to network fields, forming plages and faculae all the way down to small-scale mixed polarity turbulent magnetic field filling the rest of the solar surface. State-of-the-art solar simulations reproduce all of these features, from sunspots (Rempel et al. 2009) and plages (Vögler et al. 2005; Yadav et al. 2021) to quiet-sun magnetism (Rempel 2014).

The quiet-sun small-scale field, that is, the field associated with regions of the sun not showing any activity, is partly attributable to a small-scale dynamo (SSD) operating in the convection zone (Vögler & Schüssler 2007; Rempel 2014; Pietarila Graham et al. 2010). This field is believed to have a significant magnitude of around ~130 G, based on Hanle depolarization (Trujillo Bueno et al. 2004). Additional evidence that a fair fraction of the small-scale field is a result of an SSD comes from the fact that internetwork magnetic flux does not follow the solar cycle (Buehler et al. 2013; Lites et al. 2014). In realistic radiative-MHD simulations, the effect of quiet-sun magnetic fields (self consistently generated via an SSD mechanism) has previously been shown to be important, for example, to reproduce the correct solar intensity contrast (Danilovic et al. 2010) and account for inferred photospheric magnetic field strength based on Hanle-effect diagnostics (Shchukina & Trujillo Bueno 2011). In addition, there is a tendency to obtain a rough equipartition between kinetic and magnetic energy in SSD simulations (Hotta et al. 2015; Haugen et al. 2004; Schekochihin et al. 2004), which implies a substantial reduction in plasma velocities since most of the energy in the magnetic fields is obtained from the plasma motions. The importance of magnetic fields generated from an SSD in other stellar types and its effect on the intensity characteristics, however, remains to be explored. Hence, it is imperative to investigate the effect of quiet-star small-scale magnetic fields on these quantities in a subsequent study.

This paper is a part of a project aimed at constructing a grid of magneto-convective stellar atmospheres, ranging across temperatures (3500 < T < 7000 K), gravity (4.3 < log₁₀ g < 4.8 in cgs units), and metallicities (in this paper, only solar metallicities are considered). We investigate four stellar cases: F3V, G2V, K0V, and M0V and we study the relative change in stratification, convection, and intensity from purely hydrodynamic setups.

In Sect. 2, we outline the simulation code and the setup. Then, we present the results of the simulations in Sect. 3, followed by interpretation of the results in Sect. 4. Lastly, we summarize the results and present the corresponding discussion in Sect. 5.

2 Methods

2.1 Simulation Code

The code we use throughout this work is MURaM (Vögler et al. 2005; Rempel 2014, 2017), a 3D radiative-MHD code that solves the conservative MHD equations for compressible, partially ionized plasma. It uses a multi-group radiative transfer scheme with short characteristics (Nordlund 1982). The equations for mass (ρ − density), momentum (v − velocity, p − pressure), and energy (ϵ_h − enthalpy density) conservation are solved, along with the induction equation (B − magnetic field) $${\partial }_t\rho =-\nabla \cdot \left(\rho v\right),$$(1) $${\partial }_t\left(\rho v\right)=-\nabla \cdot \left(\rho vv\right)-\nabla p+\rho g+F_{\text{SR}}+F_L,$$(2) $$\begin{array}{r}\hfill {\partial }_t\left({\in }_h+\rho v^2/2\right)=-\nabla \cdot \left(v\left({\in }_h+p+\rho v^2/2\right)\right)\\ \hfill +v\cdot \left(g+F_L+F_{\text{SR}}\right)+Q_{\text{rad}}+Q_{\text{res}},\end{array}$$(3) $${\partial }_{t}B=\nabla \times \left(v\times B\right).$$(4)

Here, F refers to forces and Q refers to source terms. The subscript SR refers to semi-relativistic “Boris correction”-related terms (Boris 1970; Gombosi et al. 2002), which are negligible for our setups, and L refers to the Lorentz force. The two Q terms in the energy equation account for the radiative heating/cooling and resistive heating (since the hydrodynamic energy is conserved, instead of the total energy). For details, we refer the reader to Rempel (2014, 2017). In this work, the grey approximation is used for solving the radiative transfer equations, where the frequency dependence of the opacity is replaced by an average value. This is an acceptable approximation for this work since we are mainly interested in the structure below and just above the surface (Vögler et al. 2004). Finally, the FreeEOS equation of state (Irwin 2012) with solar abundances (Asplund et al. 2009) is used to close the set of equations.

The effective temperature (T_eff) (related to the radiative output), the surface gravitational acceleration (g), related to the hydrostatic balance, and the metallicity (Z) together uniquely specify the spectral class of a star. The MURaM code uses a constant gravitational acceleration g, the gas pressure at bottom boundary (p_bot), and the inflow entropy at the bottom boundary (s_bot) as free parameters. The p_bot and s_bot parameters determine the height of the τ = 1 surface and the T_eff.

2.2 Setup and Parameters

We considered four stellar cases: F3V, G2V, K0V, and M0V, chosen to cover a broad range of T_eff for stellar types with convective envelopes. All boxes have the same number of grid-points (512 × 512 in the horizontal direction and 500 in the vertical direction). The scaling for the horizontal and the vertical geometric extent was done such that the number of granules in each box is roughly the same and the number of pressure scale heights below the photosphere is also similar (~7.5). For the G-star, this corresponds to 4 Mm below the surface and a horizontal extent of 9 Mm × 9 Mm.

The boundaries are periodic in the horizontal x, y direction. The top boundary (z_top) is open to outflows and closed to inflows, with vertical magnetic fields. The bottom boundary (z_bot) is symmetric¹ for mass flux (pv), entropy downflows, and magnetic fields. This magnetic field boundary condition also allows horizontal field to be advected across the bottom boundary. This mimics” the presence of magnetic field deeper in the convection zone, as considered previously on the basis of equipartition arguments (Rempel 2014; Hotta et al. 2015). The magnetic field BC may not necessarily preserve the ∇ · B = 0 constraint. However, the hyperbolic divergence cleaning approach (Dedner et al. 2002) employed in MURaM takes care of the ∇ · B errors reasonably well: (∇ · B)_rms/(‖B‖/∆z) ~ O(10⁻³); we refer to Fig. C.1 for details.

For each star, we performed simulations with magnetic fields (SSD), and purely hydrodynamic (HD) simulations. HD simulations were initialized with density and internal energy (IE) profiles generated using the 1D stellar code MESA (Paxton et al. 2019) for the F-, K-, and M-star and using the standard solar model from Christensen-Dalsgaard et al. (1996) for the G-star. These were then run for several hours in stellar time till convection became relaxed and there were no transients visible in velocity and pressure vertical slices. Then the simulation box was seeded with net zero-flux magnetic field with a negligibly small field strength (10⁻⁵ G) and run till photospheric magnetic field strength reached saturation. The results presented in the subsequent sections are averaged over a few hours of stellar time (after saturation), and over a number of snapshots, for all the eight cases (see Table 1 for further details for each setup).

Table 2 describes the detailed setup for all the simulations: for all four stellar types, it gives the height of the τ = 1 surface above the bottom boundary (and the total vertical extent), the horizontal extent, the horizontal resolution, the vertical resolution, the log constant surface gravity, the effective temperature of the SSD and HD cases, the average magnetic field magnitude at the τ = 1 iso-surface² 〈‖В‖〉_τ=1, the average unsigned vertical field at the τ = 1 iso-surface ∣B_z∣_τ=1, and the average magnetic field magnitude at the 〈τ〉 = 1 horizontal slice 〈‖В‖〉_{〈τ〉 = 1} of the SSD cases. The effective temperature is calculated by averaging the angle-averaged bolometric luminosity over time.

3 Results

All the magnetic simulations exhibit dynamo action and develop magnetic fields with energy within an order of magnitude of the kinetic energy (KE) through most of the simulation box. The change in the partition of energy influences the internal structure as well as convective velocities. Snapshots of the bolometric intensity and vertical velocity at the τ = 1 surface for the SSD and the HD setups (along with the vertical magnetic field for the SSD cases) are shown in Fig. 1. The SSD cases show distinct intergranular bright points which correspond with strong magnetic field concentrations. In the following subsections, we examine the horizontally averaged structure of the magnetic field and its effects on the stratification as well as convection for these stellar types. The analysis of the magnetic fields in the lower photosphere and their effect on the bolometric intensity and vertical velocity will be covered in the next paper in this series.

All 1D plots are averages over a number of snapshots spanning a few hours of stellar time (see Table 1 for exact numbers). The error bars are standard error (standard deviation normalized by the square root of the number of snapshots $\sigma /\sqrt{N}$) of the average 1D structure, the assumption being that over this time span, the snapshots are statistically independent.

All quantities are plotted as a function of number of pressure scale heights relative to the height where 〈τ〉 = 1, n_H = log(p_gas/p_{gas(〈τ〉 = 1)}) We note that with this definition, positive values correspond to the interior. Since the non-magnetic bottom boundary conditions are identical (p_bot and s_bot are the same) for HD and SSD runs, the deviations are calculated geometrically and plotted against the corresponding HD pressure scale axis.

This extent in terms of pressure scale heights ranges from 7.5 (bottom) to −5 (near top) for all the simulations. Since the focus of this paper is on the near surface convection zone, we excluded the portion of the box corresponding to n_H < −1 from our analysis. We also excluded the region corresponding to n_H > 6 due to possible numerical bottom boundary effects.

3.1 Magnetic Field Structure

As mentioned in the introduction to this section, the dynamo action results in a significant amount of magnetic field, with the overall magnitude roughly similar for all the cases, and with a somewhat decreasing trend with T_eff near the surface (Fig. 2, top panel and last column of Table 2). The relation of the magnetic energy (ME) to kinetic energy is discussed in Sect. 3.3.

The bottom panel of Fig. 2 shows the ratio B_h,rms/B_z,rms, which gives an indication of the 3D structure of the magnetic fields. For a fully isotropic distribution of magnetic field, one would expect $B_x^2\approx B_y^2\approx B_z^2$. This implies that the ratio of the horizontal r.m.s. component of the magnetic field fluctuations B_h,rms to the vertical component B_z,rms should be $\sim \sqrt{2}$. For G-, K-, and M-stars, this ratio is slightly less than $\sqrt{2}$ in the middle of the box (2 < n_H < 5), indicating near-isotropy, whereas for the F-star, it is significantly lower (~1).

Near the τ = 1 surface, B_h,rms/B_z,rms for all the stars is lower because of intensification of vertical magnetic fields in the intergranular lanes (e.g., Spruit 1979).

3.2 Changes in Stratification

Figure 3 shows the plots of deviation from the mean HD stratification. The deviations from the HD simulations are presented as the relative percent change in the horizontally averaged 1D structure. For any quantity of interest q (e.g., density, temperature, etc.), these deviations are calculated as (q_SSD − q_HD)/q_HD. This means that a positive value for the deviation corresponds to a higher value for the SSD case relative to the HD case.

All simulations show slight (≤2%) changes in thermodynamic stratification relative to the corresponding HD simulations, with the magnitude of deviations below the surface roughly increasing from coolest to hottest stellar type. Below 3 pressure scale heights, the G, K, and M simulations show negligible (≤0.1%) deviations in thermodynamic quantities.

Closer to the surface, the F-star simulation shows up to 1.5% reduction in density and up to 1% reduction in gas pressure (Fig. 3, top and bottom panels, blue line). This trend is opposite to that seen for other cases, all of which show a slight (<1%) enhancement in density and pressure. These results are analyzed in Sect. 4.1. For context, a 1.5% deviation in temperature for an F-star would correspond to a temperature change of ~ 100 K.

3.3 Distribution of Energies

In Fig. 4, panel a shows the ratio of magnetic to kinetic energy. For all cases, this ratio is within an order of magnitude throughout the box. This ratio for the F-star is significantly higher compared to the other stars. In addition, the ratio for M-star has a minima deeper down in the box than for other stars. As before, the general trend shows a decrease in this ratio with decreasing T_eff.

In Fig. 4b, we shows the relative change in KE between the SSD and the HD setups. All stars show a marked decrease in KE in the SSD case. As with the top panel, F-star shows a significantly stronger reduction in KE compared to the other cases. The M-star as well shows a stronger reduction in KE near the surface compared to the G and K-stars.

Figure 4c shows the ratio of kinetic to internal energy. For the F-star, this ratio is within an order of magnitude near the surface, whereas for the other stars, it is less than 1%. Figure 4d shows the ratio of the magnetic to the internal energy. There is a clear trend in this ratio with stellar type, with it being lowest for M-star (~5 × 10⁻⁵ near bottom, to ~10⁻⁴ near surface) and highest for F-star (~10⁻³ near bottom, to ~10⁻² near surface).

Fig. 1 Emergent intensity and surface vertical velocities in different stellar types for models with and without magnetic field. From top to bottom: snapshot of the bolometric intensity and vz at τ = 1 for the HD case (rows 1 and 2), bolometric intensity and vz at τ = 1 for the SSD case (row 3 and 4), and the corresponding vertical magnetic field at τ = 1 (row 5, from left to right) for spectral types F, G, K, and M, respectively. The green circles indicate the bright points and corresponding magnetic field concentrations.

Fig. 2 Magnetic field structure for the four stellar cases. Top: horizontally averaged magnetic field magnitude. Bottom: the ratio of the horizontal rms field strength to the vertical r.m.s field strength. The horizontal axis is the number of pressure scale heights log10(pgas/pgas(τ=1)), calculated for the HD cases, below the surface (dotted vertical black line). The shaded regions correspond to 1-σ standard error N is the number of snapshots) computed over time averaging of snapshots.

Fig. 3 Deviations in ρ, T, and pgas for F (blue), G (black), K (green), and M-star (red) cases. The vertical axis gives the geometric deviations as a percentage relative to the hydrodynamic case. The horizontal axis is the number of pressure scale heights log10(pgas/pgas(τ=1)), calculated for the HD cases, below the surface (dotted vertical black line). The shaded regions correspond to 1-σ standard error $(e=\sigma /\sqrt{N}$, N is the number of snapshots) of the mean solid curve.

Fig. 4 Distribution of energies for the four stellar cases. (a) Ratio of ME to KE, (b) percent change in KE, (c) ratio of KE to IE, and (d) ratio of ME to IE (bottom) for F-, G-, K-, and M-stars.

3.4 Changes in Velocities

All SSD cases exhibit a decrease in vertical velocities $v_z^2$ as well as the ratio of horizontal to vertical rms velocities v_h,rms/v_z,rms (Fig. 5), relative to the corresponding HD cases. The reduction in $v_z^2$ follows similar trend for all the four stars, with a decrease of 4–8% near the surface, going up to 20% near the bottom boundary. For G-, K- and M-star, the decrease in v_h,rms/v_z,rms is similar (5% near the surface, going up to 12% near the bottom) but is more pronounced for the F-star case (10% near the surface, going up to 20% near the bottom boundary). This implies a change in the horizontal extent of sub-surface granulation, which follows Nordlund et al. (2009), who showed using simple mass conservation that the horizontal extent of granules is proportional to H(v_h/v_z), where H is the local density scale height.

The ratio v_h,rms/v_z,rms gives an idea of the 3D velocity structure. Close to the τ = 1 surface, the ratio increases suddenly. This increase corresponds to where convective flows turn over, as the atmosphere becomes convectively stable and the flows above that point are mainly due to convective overshoot. The exact pressure scale depth where this turning over takes place depends of the effective temperature: for the F- and G-star, this turnover takes place within half a pressure scale height of the τ = 1 surface whereas, for the M-star, it takes place well below the τ = 1 surface (around n_H = 3). This is also probably why the minima of the ME/KE ratio in Fig. 4 for M-star is significantly below the surface compared to the G and K-star, as density (and, consequently, KE) is lower above the height where most of the overturning takes place. The changes in $v_z^2$ as well as the v_h,rms/v_z,rms with depth seem to follow a regular trend, and is consistent with the results obtained previously in whole-convection zone simulations with a small-scale dynamo setup (Hotta et al. 2015).

Fig. 5 Relative decrease in convective velocities $v_z^2$ (top), relative decrease in ratio of horizontal to vertical rms velocities vh,rms/vz,rms (middle) and the actual vh,rms/vz,rms ratio (bottom).

4 Discussion

4.1 Turbulent Pressure

We first discuss the changes in the pressure and density in the SSD models relative to a purely HD model. The pressure and density changes (Fig. 3) are very small (~0.1%) for the G-, K-, and M-star in the convection zone (and not significantly greater near the surface for the K- and M-star). The changes are more prominent for the F-star This discrepancy can be understood by considering the contribution of turbulent pressure in the overall hydrostatic balance.

Turbulent pressure becomes important for hydrostatic balance when KE is within an order of magnitude or so of the IE. This is indeed the case for our F-star models. This implies that turbulent pressure (which is largely due to plasma motions) can be a significant fraction of the gas pressure p_gas. Quantitatively, this can be seen from a crude MLT calculation of the Mach number (see Appendix A for derivation): $$M\approx 0.138{\left(\frac{T}{{10}^3\text{K}}\right)}^{5/6}{\left(\frac{\rho }{{10}^{-7}{\text{g cm}}^{\text{-3}}}\right)}^{1/3}.$$(5)

Using the above equation, we calculate the F-star photospheric Mach number to be about 0.75 (see Table 3), slightly lower than simulation value in Beeck et al. (2013) who found it to be 0.9. For the other stars, the velocity is decidedly subsonic.

Now, the hydrostatic balance is expressed in terms of the balance between force due to pressure gradient, p′ = dp/dz, and gravity, ρg, acting in the vertical (viz. radial) direction. These two terms should be approximately equal in magnitude. If just p_gas is considered, this balance does not hold very well, with deviations increasing strongly with T_eff (Fig. 6, top panel).

The turbulent pressure consists of terms from the total stress tensor. From Reynolds stresses, the ρv_iv_j term and from Maxwell stresses, the ((B²/2)δ_ij − B_iB_j)/4π term is obtained. As mentioned in the introduction, the presence of SSD magnetic fields implies a reduction in KE, where energy is redistributed between the plasma motions and magnetic fields via Lorentz force feedback. In Fig. 4, the ratio of ME to KE (panel a) corresponds remarkably well with the reduction in KE relative to the HD case (panel b), implying most of the energy in the magnetic fields is drawn from the KE reservoir. The ME is within an order of magnitude of the KE for subsurface plasma. This is consistent with the results on equipartition of energy in several SSD simulations (Hotta et al. 2015; Haugen et al. 2004; Schekochihin et al. 2004).

This order of magnitude equipartition results in a reduction of plasma velocities and, consequently, the magnitude of the Reynolds stress term. The contribution from the Maxwell stress, on the other hand, can be either negative or positive, depending on whether the effects of magnetic tension dominate over that of magnetic pressure in the vertical direction. More quantitatively (see Appendix B for a derivation), the total turbulent pressure can be expressed as: $$p_{\text{turb}}=\rho v_z^2+\frac{B_{\text{h}}^2-B_z^2}{8\pi }.$$(6)

With this included in the calculation of the pressure gradient term, the hydrostatic balance is satisfied, as can be seen in middle panel of Fig. 6 (noting the difference in the extent of y-axis for the top and middle panels). In the bottom panel of the same figure, the reduction in p_turb /p_gas is most prominent for the F-star. Since the changes in $v_z^2$ are similar (Fig. 5, top panel) for all stars, there must be a reduction in density for the F-star to compensate for the significant change in p_turb. Based on the expression for p_turb, we introduce an effective turbulent velocity $v=\sqrt{p_{\text{turb}}/\rho }$.

With this, it becomes possible to relate changes in density stratification to changes in total pressure gradient. A straightforward consideration of the density scale height H_ρ (see Appendix B) yields: $$H_{\rho }=RT/\left(\mu g\right)+v^2/g.$$(7)

This tells us that where $B_{\text{h}}^2<B_{\text{z}}^2$ and ${\left(\rho v_z^2\right)}_{\text{SSD}}\text{<}{\left(\rho v_z^2\right)}_{\text{HD}}$, the scale height of the SSD model is smaller than for the HD model, (H_ρ)_SSD < (H_ρ)_HD. This is valid if assuming that the change in T and μ is relatively small. From the proportionality between density scale height and $v_z^2$ (Eq. (7)), the decrease in vertical velocity is associated with a decrease in the local density scale height. This can be inferred from Fig. 6 (bottom panel), where the ratio p_turb/p_gas (∝ v²/(RT/μ)) is noticeably lower for F-star near the surface.

These conditions are satisfied relatively well for the F-star case. A decrease in H_ρ implies a steepening of the density stratification for the SSD case relative to the HD case. This is exactly the case in the top panel of Fig. 3: the reduction in density goes from a ~0% to ~1.5% near the surface. Below the surface, it is reasonably good to assume a perfect gas equation of state, with p_gas ∝ ρT. Assuming changes in T are small compared to ρ and p_gas, this, in turn, implies that the changes in p_gas follows a similar trend as changes in ρ. Near the surface, this picture breaks down since the mode of energy transfer changes from convective to radiative and this simplified analysis is no longer valid.

The effect of turbulent pressure on stratification and convection has previously been considered for stellar model envelopes from an MLT perspective (Henyey et al. 1965) as well as for 3D HD simulations (Ludwig & Kučinskas 2012; Jørgensen & Weiss 2019); however, to the best of our knowledge, the effect of magnetic fields from an SSD have not been considered before.

Fig. 6 Hydrostatic force balance with gas pressure pgas (top) and total pressure ptot (middle) gradient term and ρg term, normalized by 〈ρg〉 as 〈dp/dz + ρg〉/〈ρg〉, and the ratio of turbulent pressure to gas pressure (bottom). All quantities are plotted for SSD (solid) and HD (dashed) cases. We note that the extent of y-axis is larger by more than an order of magnitude for force balance with pgas compared to force balance with ptot.

4.2 Changes in Velocity Structure

As mentioned in Sect. 3.4, the ratio v_h,rms/v_z,rms (Fig. 5, bottom panel), gives an idea of where convection overturns relative to pressure scale height for a given effective temperature. The trend in overturning for different spectral types fits the usual picture of “hidden” granulation below the τ = 1 surface for cooler, denser stars such as M and K versus the “naked” granulation above the τ = 1 surface of hotter, more rarefied stars such as G and F, as discussed in Nordlund & Dravins (1990).

This ratio is close to 1 for the HD G-, K-, and M-stars (dashed lines). For the HD F-star, however, it is significantly higher (dashed light blue line), indicating a higher degree of isotropy in the velocity structure³. This implies relatively higher horizontal velocities for the HD F-star. The near-isotropic velocity profile for the HD F-star can be attributed to the KE being a non-negligible fraction of the IE. We speculate that these stronger horizontal velocities contribute to the stronger magnetic fields in the vertical direction for the F-star (as inferred from the lower B_h,rms/B_z,rms for the F-star in Fig. 2, bottom panel).

In the case of SSD (solid lines), all stars follow a similar trend for v_h,rms/v_z,rms, with the ratio being ≤1 up to the point where convection overturns. The change in velocity field structure in the presence of SSD fields is a hard problem that depends on the SSD saturation mechanism. We note that the B_h,rms/B_z,rms ratio is lower for F-star compared to other cases. It is plausible that the relatively stronger fields in the vertical direction restrict horizontal flows more in the case of F-stars and lead to a greater change in the value of v_h,rms/v_z,rms. However, a fuller understanding of this behavior requires a more detailed analysis that is beyond the scope of this paper.

5 Conclusions

In this work, we investigate the magnetic field self-consistently generated by an SSD acting in the near-surface layers of main-sequence stars of spectral types F3V, G2V, K0V, and M0V. The SSD mechanism operates in all cases to amplify magnetic fields from a seed field of negligible strength and zero net flux. The magnetic fields from the SSD have an energy density that is a non-negligible fraction of the kinetic energy density. These fields act back on the plasma to reduce the convective velocities, which in turn reduces the turbulent pressure. This becomes substantial for the F-star as it is hot enough to have kinetic and internal energy within an order of magnitude near the surface, which gives magnetic fields stronger than those in G-, K-, and M-stars, especially in the vertical direction. The equation for hydrostatic balance for total pressure and the reduction of convective velocities implies a reduction in the density scale height itself. This is significant enough for the F-star to result in reduced density and gas pressure throughout the box. This effect tends to get smaller towards later spectral types.

This paper only covers the near-surface convection zone. Other aspects of particular interest to observational studies would include the magnetic field structure in the lower photosphere and changes in the intensity characteristics. In addition, the change in scale height and the changes in v_h/v_z imply changes in granulation scale. All these points are to be investigated in a follow-up paper.

Acknowledgements

T.B. ran the simulations, performed the analysis and wrote most of the text. D.P. helped with running the simulations. All authors contributed to the final version of the text. We thank the anonymous referee for their helpful and encouraging comments. T.B. is grateful for access to the supercomputer Cobra at Max Planck Computing and Data Facility (MPCDF), on which all the simulations were carried out. We also acknowledge Alan Irwin for providing the open source FreeEOS equation of state. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement nos. 695075 and 715947).

Appendix A Mach number derivation

We derive the expression for Mach number used in Eq. (5) from §4.1. Considering the hydrostatic pressure balance equation, the expression for pressure-scale height (H_p), the expression for sound speed (c_s) and the ideal gas equation of state: $$p\text{'}=-g\rho ,$$(A.1) $$H_p=-p/p\text{'},$$(A.2) $$c_s^2=\gamma p/\rho ,$$(A.3) $$p=\rho \left(R/\mu \right)T.$$(A.4)

Here, H_p is the pressure scale height, p′ = dp/dz, μ is the mean molar mass, and 1/μ = (1 + E)(X + Y/4 + Z/2), where E is ionization fraction and X, Y, Z are the H, He and metal abundances. In the Böhm-Vitense MLT (Böhm-Vitense 1958), convective flux and velocities are expressed as: $$F_{conv}=\alpha c_p\rho vT\left(\nabla -{\nabla }_a\right)/2,$$(A.5) $$v^2={\alpha }^2\delta g{H}_p\left(\nabla -{\nabla }_a\right)/8.$$(A.6)

Here, α is the mixing length parameter (usually taken to be between 1.5 to 2), δ = 1 − (∂ ln μ/∂ ln T)_p (abundance gradient with temperature for an ideal gas), c_p = (∂U/∂T)_p + pδ/(ρT) (heat capacity at constant pressure), and ∇ = d ln T/d ln p, with the subscript a referring to ‘adiabatic.’ We refer to chapter 6 of Stix (2004) for more details.

Near the surface, as most energy is carried by radiation, F_conv ≈ σT⁴. We take the ratio of Eq. A.5 and Eq. A.6 to eliminate ∇ − ∇_a and use σT⁴ instead of F_conv. Taking Mach number M = v/c_s and $c_s=\sqrt{\gamma RT/\mu }$ and eliminating H_p and p using Eq. A.1, A.2, and A.4, $$M={\left(\frac{\alpha \delta \sigma {\mu }^{1/2}}{4c_p{\gamma }^{3/2}{R}^{1/2}}\right)}^{1/3}{T}^{5/6}{\rho }^{-1/3}.$$(A.7)

Using relevantvalues forthe terms intheparentheses (allincgs), we get Eq (5) in terms of T and ρ as above:

The value of μ = 1.2 assumes E ~ 0 near τ = 1. For an M-star, this may not necessarily be the case. However, even if E → 1, μ changes by a factor of 2 and, accordingly, M changes by a factor of 2^1/2 ≈ 1.4, and the qualitative result still holds.

Appendix B Turbulent pressure and pressure scale height derivation

In §4.1, the expression used for p_turb in Eq. (6) can be derived by considering the total (Reynolds and Maxwell) stress tensor, σ_ij, for ideal MHD: $${\sigma }_{ij}\rho \left(v_i{v}_j\right)+\left(p+\frac{B^2}{8\pi }\right){\delta }_{ij}-\frac{B_i{B}_j}{4\pi }.$$(B.1)

Here i, j represent the x, y and z directions.

The pressure balance along the vertical (z) direction involves the term 〈∇_Z · σ_iz〉_z which is equal to ∂p_tot/∂z. We can now compare these two terms to get an expression for total pressure (p_tot). Assuming all the off-diagonal terms (i ≠ j) are negligible (which is akin to saying there is no cross-correlation between the vertical and the horizontal velocity and magnetic field components), we obtain the following: $${\langle p_{\text{tot}}\rangle }_z=\langle p_z^2\rangle +{\langle p\rangle }_z+\frac{\langle B_h^2-B_z^2\rangle }{8\pi }.$$(B.2) $${\langle p_{\text{tot}}\rangle }_z=\langle p_z^2\rangle +{\langle p\rangle }_z+\frac{\langle B_h^2-B_z^2\rangle }{8\pi }.$$(B.3)

The next step is to consider how using p_tot instead of p affects the density scale height, H_ρ. In Eq. (A.1), we substitute expression for p_tot, instead of p, and use Eq. (A.4) to eliminate p. With this, we obtain: $$\frac{d}{dz}\left(\rho \left(RT/\mu +v^2\right)\right)=-\rho g.$$(B.4)

Here, we take v² to be $v_z^2+\left(B_h^2-B_h^2\right)/\left(8\pi \rho \right)$. Then, the above equation can be rearranged to obtain: $$\frac{d\rho }{\rho }=-\frac{gdz}{RT/\mu +v^2}-\frac{d\left(RT/\mu +v^2\right)}{RT/\mu +v^2}.$$(B.5)

From this, we can obtain the expression for density scale height, H_ρ, as: $$H_{\rho }=RT/\left(\mu g\right)+v^2/g.$$(B.6)

The second term above includes the contribution from the turbulent pressure. Since this term is almost always smaller for the SSD case relative to the HD case (mainly because the convective velocities are lower for SSD), the corresponding pressure scale height for the SSD cases is also almost always smaller.

Appendix C Diagnostics

For any MHD simulation, the computation and evolution of the magnetic field must be divergence free. To ensure this, MURaM uses a hyperbolic divB cleaning algorithm (Dedner et al. 2002). Here, we show the horizontally averaged divergence of magnetic field across the box for all four SSD cases. Since the units for ∇ · B are field/length, a proper comparison requires a normalization. We do so with ‖B‖/∆z. As noted in §2.2, at best, the error is O(10⁻³).

Fig. C.1 div · B error for the four SSD cases: F, G, K, and M. The quantity plotted is horizontally averaged (∇ • B)rms//(‖B‖/∆z..

References

All Tables

All Figures

Fig. 1 Emergent intensity and surface vertical velocities in different stellar types for models with and without magnetic field. From top to bottom: snapshot of the bolometric intensity and vz at τ = 1 for the HD case (rows 1 and 2), bolometric intensity and vz at τ = 1 for the SSD case (row 3 and 4), and the corresponding vertical magnetic field at τ = 1 (row 5, from left to right) for spectral types F, G, K, and M, respectively. The green circles indicate the bright points and corresponding magnetic field concentrations.

In the text

Fig. 2 Magnetic field structure for the four stellar cases. Top: horizontally averaged magnetic field magnitude. Bottom: the ratio of the horizontal rms field strength to the vertical r.m.s field strength. The horizontal axis is the number of pressure scale heights log10(pgas/pgas(τ=1)), calculated for the HD cases, below the surface (dotted vertical black line). The shaded regions correspond to 1-σ standard error N is the number of snapshots) computed over time averaging of snapshots.

In the text

Fig. 3 Deviations in ρ, T, and pgas for F (blue), G (black), K (green), and M-star (red) cases. The vertical axis gives the geometric deviations as a percentage relative to the hydrodynamic case. The horizontal axis is the number of pressure scale heights log10(pgas/pgas(τ=1)), calculated for the HD cases, below the surface (dotted vertical black line). The shaded regions correspond to 1-σ standard error $(e=\sigma /\sqrt{N}$, N is the number of snapshots) of the mean solid curve.

In the text

	Fig. 4 Distribution of energies for the four stellar cases. (a) Ratio of ME to KE, (b) percent change in KE, (c) ratio of KE to IE, and (d) ratio of ME to IE (bottom) for F-, G-, K-, and M-stars.
In the text

	Fig. 5 Relative decrease in convective velocities $v_z^2$ (top), relative decrease in ratio of horizontal to vertical rms velocities vh,rms/vz,rms (middle) and the actual vh,rms/vz,rms ratio (bottom).
In the text

Fig. 6 Hydrostatic force balance with gas pressure pgas (top) and total pressure ptot (middle) gradient term and ρg term, normalized by 〈ρg〉 as 〈dp/dz + ρg〉/〈ρg〉, and the ratio of turbulent pressure to gas pressure (bottom). All quantities are plotted for SSD (solid) and HD (dashed) cases. We note that the extent of y-axis is larger by more than an order of magnitude for force balance with pgas compared to force balance with ptot.

In the text

	Fig. C.1 div · B error for the four SSD cases: F, G, K, and M. The quantity plotted is horizontally averaged (∇ • B)rms//(‖B‖/∆z..
In the text